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propagation”  network,  and  a  network  which  learns  by  a  method  of 
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networks  as  cognitive  systems  and  neural  networks  as  mere  associ¬ 
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BIOMASSCOMP  PHASE  I  FINAL  REPORT 


1.  INTRODUCTION  AND  OBJECTIVES 

This  is  a  technical  report  of  a  six  month  Phase  I  research 
project  supported  by  the  U.S.  Air  Force  Wright  Aeronautical 
Laboratories  {WPAFB,  OH)  under  the  Small  Business  Innovation 
Research  Program.  This  report  presents  a  complete  account  of  our 
investigations  of  the  subject  pursuant  to  the  objectives  of  the 
original  SBIR  proposal,  and  is  organized  according  to  the  listing 
of  those  objectives  in  the  proposal.  In  the  pursuit  of  those 
objectives,  we  have  not  only  achieved  the  implementation  of  a 
successful  entropic  index  of  the  structural  similarity  of  two 
systems,  but  we  have  also  identified  errors  in  our  planned 
approach  through  the  conduct  of  experiments  that  failed  to  con¬ 
clusively  demonstrate  the  expected  results.  This  report  details 
both  the  successes  and  the  failures,  and  the  valuable  information 
that  we  have  learned  from  them. 

As  stated  in  our  Phase  I  proposal,  our  objective  was  to 
demonstrate  the  feasibility  of  developing  and  applying  an 
entropic  measure  of  structural  similarity  of  systems  so  as  to 
obtain  (in  the  follow-on  project)  an  automated  procedure  for 
mapping  the  architecture  of  a  living  neural  network  into  a 
machine.  In  order  to  accomplish  this  objective,  our  tasks  were; 

1.  Identify  and  develop  a  mathematical  technique  for  the 
measurement  amd  analysis  of  relative  information  content 
in  the  signals  of  a  network, 

2.  Identify  and  develop  a  mathematical  technique  for  the 
parametric  optimization  of  artificial  neural  network 
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models  as  measured  by  the  combined  relative  information 
content  of  a  functioning  hybrid  network, 

3.  Identify  the  functional  design  of  a  multichannel  bidir¬ 
ectional  signal  translator  suitable  for  the  realtime 
interface  of  a  natural  network  on  the  multimicioelec- 
trode  plate  (MMEP)  apparatus  of  G.  Gross  at  North  Texas 
State  University  (NTSU)  to  an  artificial  network, 

4.  Closely  monitor  euid  assist  the  ongoing  work  at  NTSU  to 
demonstrate  the  capability  of  extracellular  electrodes 
in  the  MMEP  apparatus  to  be  used  for  the  injection  of 
localized  potentials  capable  of  stimulating  activity  in 
specific  subnets  of  the  cultured  neural  network, 

5.  Closely  monitor  and  assist  the  ongoing  work  at  NTSU  to 
demonstrate  the  ability  to  "condition"  the  behavior  of  a 
natural  neural  network  in  culture  through  controlled 
stimulation. 

Our  principal  achievement  has  been  the  definition  and 
algorithmic  implementation  of  a  scalar  measure  of  structural  sim¬ 
ilarity  of  two  systems,  based  on  extensive  time-series  of  state 
measurements  (signal  vectors)  from  the  two  systems.  This  report 
details  that  definition,  and  the  FORTRAN  source  code  of  the 
algorithm  is  included  in  an  appendix.  A  series  of  experiments 
with  random  data  vectors  containing  varying  degrees  of 
correlations  demonstrates  the  behavior  of  the  algorithm. 

Moreover,  these  experiments  predict  an  interesting  psychological 
tendency  of  an  observer  to  always  perceive  the  greatest  degree  of 
order  in  his  observed  environment  when  his  knowledge  is  at  a 
certain  intermediate  stage  between  total  ignorance  and  complete 
understanding. 

The  application  of  DMORPH  to  neural  network  architectures 
has  shown  that  our  approach  to  using  the  structure  measure  to 
improve  the  architectural  design  of  neural  networks  by  comparing 
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them  natural  networks  appears  to  have  been  flawed.  Yet,  the 
flaw  is  not  one  that  might  have  been  easily  detected  without  a 
study  of  the  results  of  the  experiments  (although  the  a  posteri¬ 
ori  explanation  in  the  form  of  a  suitable  "gedankenexperiment"  is 
simple  enough),  and  its  exposure  has  resulted  in  new  emd  useful 
insights  into  the  structural  and  functional  principles  of  neural 
networks  and  other  self-organizing  systems.  Those  insights  are 
discussed  in  the  "Analysis"  section  of  this  report,  and  they  will 
constitute  the  direction  for  our  planned  Phase  II  research. 

1 . 1  Background 

The  outline  of  our  argument  is  this:  The  objective  of  an 
intelligent  system  is  to  minimize  surprise,  or  novelty,  in  its 
interaction  with  its  environment  in  a  manner  that  is  consistent 
with  its  "mission".  To  this  end,  it  builds  predictive  models  of 
the  world  and  stores  these  models  in  any  convenient  recording 
medium  including,  but  not  limited  to,  its  own  memory.  A  predic¬ 
tive  model  will  be  more  carefully  defined  below  in  Section  7,  but 
heuristical ly  it  is  a  transition  operator  which  associates  each 
sensory  measurement  with  an  empirically-based  probability  density 
for  the  perceptive  effects  of  future  observations.  (Ho  &  Lee 
(4])  This  description  is  further  illustrated  by  Watanabe,  who 
says  (CIO],  p.l42)  "The  existence  of  structure  means  that  the 
knowledge  of  a  part  allows  us  to  guess  easily  the  rest  of  the 
whole. 

The  brains  of  humans  and  animals  are  not  apart  from  the 
universe,  but  are  parts  of  the  whole.  Their  function,  which  we 
summarize  with  the  verb,  "to  learn",  is  to  adopt  a  form  which, 
when  explored  by  the  animal  through  associative  recall,  allows  it 
to  guess  what  is  going  on  in  the  rest  of  the  universe  and  to 
adjust  its  behavior  to  minimize  surprize  subject  to  its  mission. 
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The  design  of  artificial  neural  networks  and  neurocomputers 
normally  proceeds  by  using  the  best  available  data  from  the 
neuroscience  community  to  build  and  test  computational  models  of 
the  components  and  structures  of  the  brain.  Models  that  work 
well  are  improved  upon.  Models  that  flop  are  filed  under 
"experience".  The  BIOMASSCOMP  project  (originally  described  in 
Dawes  [2])  was  designed  to  speed  up  this  developmental  process  by 
defining  a  computable,  numerical  measure  of  the  quality  of  a 
neural  network  model.  This  measure  estimates  the  degree  to  which 
two  neural  networks  are  producing  signals  that  have  the  same 
structure.  A  low  value  of  the  measure  means  there  is  little 
similarity  in  structure  between  the  two  systems.  A  high  value 
means  the  structures  are  strongly  correlated. 

The  initial  application  of  the  structure  function  was 
visualized  to  work  as  follows:  An  artificial  and  a  natural 
neural  network  would  be  connected  as  in  Figure  1  by  a  bidirec¬ 
tional  communication  link  called  the  "synthetic  axon  bundle". 

The  signals  emanating  from  the  natural  network  would  be  demodu¬ 
lated  with  a  pulse-rats  demodulator  and  would  constitute  one  of 
two  signal  vectors.  The  signals  emanating  from  the  artificial 
network  would  constitute  the  other  signal  vector.  The  structure 
function  would  be  applied  to  these  two  signal  vectors  and  a 
numerical  "structural  similarity"  would  be  obtained.  As  the  two 
networks  adapted  through  their  learning  laws,  we  would  expect  to 
see  changes  in  the  similarity  value  due  to  the  intrinsic  self- 
organizational  behavior  of  both  networks,  and  we  would  expect  to 
see  this  value  stabilize  asymptotically  in  the  absence  of 
external  stimuli  to  either  system.  Then,  if  we  made  any 
adjustment  to  the  architectural  parameters  of  either  system,  we 
would  expect  to  see  the  structure  value  change  again  and 
stabilize  on  a  different  value.  We  would  infer  that  the  higher 
value  of  the  structure  function  was  obtained  with  the  better  set 
of  architectural  parameters  and  we  could  therefore  obtain  further 
improvements  in  the  architecture  by  adjusting  the  parameters  in 
the  direction  of  the  highest  structural  similarity.  Since  the 
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design  parameters  of  the  simulated  network  are  under  the  dynamic 
control  of  the  experimenter’s  computer  program,  it  can  automatic¬ 
ally  increment  or  decrement  these  parameters  so  as  to  drive  the 
value  of  the  structure  function  to  its  highest  value. 


Figure  1.  The  Hybrid  Artificial /Natural  Neural  Network 


Aside  from  the  problems  illuminated  later  by  our  experi¬ 
ments,  there  are  a  number  of  difficulties  that  we  could  and  did 
anticipate.  One  of  these  is  that  the  structure  function  is 
difficult  to  compute,  but  we  have  made  considerable  progress  in 
that  respect,  as  we  shall  demonstrate.  More  serious  is  the  fact 
that  the  performance  of  complex,  nonlinear  systems  does  not 
always  improve  or  degrade  continuously  as  a  function  of  their 
design  parameters.  In  particular,  the  performance  may  be  subject 
to  bifurcations,  catastrophes,  and  chaotic  behavior  as  certain 
parameters  approach  critical  values.  This  means  that  the  search 
for  improvement  may  have  to  be  undertaken  through  stochastic 
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methods,  such  as  simulated  annealing,  rather  than  by  the  more 
standard  "hill-climbing”  methods.  Not  the  least  of  the  problems 
is  the  realization  of  the  bidirectional  communication  link 
between  an  artificial  and  a  natural  neural  network,  which  relies 
on  the  successful  development  of  a  method  for  multiple-site 
stimulation  of  the  natural  network  in  the  MMEP  culture. 
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2.  DEFINITION  OF  THE  STRUCTURE  MEASURE 

Our  initial  task  has  been  to  develop  the  measure  of  struc¬ 
tural  similarity  of  the  two  networks,  based  on  the  concept  of  the 
Gibbs  relative  entropy  function  as  described  by  Watanabe  [101. 

We  now  have  a  structure  estimator  {called  DMORPH)  running  and 
have  tested  it  on  actual  pre-processed  data  from  Prof.  Guenter 
Gross’s  laboratory  at  North  Texas  State  University  (cf..  Appendix 
A)  . 


Any  neural  network  which  is  worth  its  salt  does  at  least 
this  one  job  well:  It  builds  internal  representations  of 
external  events  in  such  a  v/ay  that  an  appropriate  stimulus  at  a 
later  time  will  recover  "a  substantial  portion"  of  the  entire 
representation.  Some  have  referred  to  this  behavior  as  "self¬ 
organization",  but  that  is  a  seriously  misleading  phrase. 

Consider  the  following  homely  analogy:  Two  politicians,  Mr.  A 
and  Mr.  B,  have  widely  differing  world-views.  Mr.  A  sees 
everything  as  either  black  or  white,  with  nothing  in  between. 

Mr.  B  maintains  a  complex  set  of  concepts  and  classifications  for 
analyzing  events.  A  third  politician,  Mr.  Z,  has  just  died. 

In  terms  of  organization,  as  expressed  by  the  entropy  level 
of  his  neural  activation  states,  Mr.  Z  is  extremely  well- 
organized,  since  the  activation  state  of  his  neurons  is  a  delta 
function  in  time  and  space.  (This  is  to  be  distinguished  from 
his  thermod^Tiamic  entropy  which,  while  lower  than  that  of  a 
warmer  living  brain,  is  still  rather  high. )  Mr.  A  is  almost  as 
well  organized,  his  activation  state  falling  always  into  one  of 
two  event  categories,  regardless  of  the  facts  of  the  external 
world.  But  Mr.  B  has  the  entropy  of  a  very  disorganized  person, 
since  his  mental  state  can  be  found  in  any  of  a  large  number  of 
configurations  over  time.  Unfortunately,  if  he  is  a  member  of 
the  "v,'rong"  political  party,  Mr.  B's  ideas  may  reflect  no  corre¬ 
lation  whatsoever  with  reality  in  spite  of  their  complexity. 
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Thus  we  see  that  organization,  as  measured  by  entropy,  is 
not  to  be  mistaken  as  an  objective  function  for  intelligence. 
There  is  no  immediate  monotonic  relationship  between  the  entropy 
(of  the  probability  distribution)  of  the  state  of  one's  thoughts 
and  the  elusive  quality  we  call  intelligence.  Yet  the  concept  of 
entropy,  properly  applied,  can  help  us  to  make  this  quality  much 
less  elusive. 

The  following  sections  detail  the  background  and  development 
of  DMORPH.  Additional  technical  details  can  be  found  in  Jaynes 
[5]  and  in  Watanabe  [10]. 


2.1  Introduction  to  Entropy. 

Entropy  is  a  real-valued  function  whose  domain  is  the  set  of 
probability  measures  on  some  given  probability  space.  When  a 
probability  measure  has  a  Radon-Nikodym  derivative,  p(x),  this  is 
called  the  probability  density  function  (pdf)  of  the  probability 
measure,  and  in  this  case,  the  entropy  of  p(x)  is  defined  by 


E(p)  =  p(x) log[p(x) 3dx  . 

Whenever  the  probability  measure  is  finite,  i.e.,  there  are  only 
a  finite  number  of  events  covering  the  sample  space  of  x,  then 
the  integral  above  can  be  replaced  by  the  sum 

E(p)  =  -  SUM!  p^  log(p^)  }  (1) 

over  the  event  sets  with  nonzero  probabilities,  p^  (which  we 
shall  refer  to  as  the  "nontrivial  events  of  p").  It  is  not 
difficult  to  see  that  this  quantity  can  range  between  a  minimum 
value  of  zero,  and  a  maximum  value  of  log(N),  where  N  is  the 
number  of  distinct  nontrivial  events  of  p.  The  minimum  value  is 
taken  when  there  is  one  certain  event  (in  which  case  N=l).  The 
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maximum  value  is  taken  when  =  1/N  for  each  nontrivial  event 

of  p. 

In  the  following,  we  shall  often  speak  of  "the  entropy  of  a 
system".  This  will  always  mean  the  entropy  of  the  probability 
density  function  for  the  states  that  the  system  can  take.  Of 
course,  the  set  of  states  for  a  given  system  is  itself  an 
abstraction  and  one  may  use  different  state  spaces  for  different 
purposes.  For  our  purposes,  we  are  interested  in  the  activation 
state  of  an  ensemble  of  neurons,  and  not  in  their  molecular 
kinetics. 

The  entropy  of  a  neural  system  is  an  extensive  quantity. 

That  is,  if  the  system  is  partitioned  into  two  subsystems,  there 
will  be  three  possibly  distinct  entropies  to  deal  with:  Those  of 
the  two  subsystems,  and  that  of  the  whole  system.  Watanabe  shows 
us  how  to  relate  these  three  quantities. 


2.2  Entropic  Structure 

Suppose  that  we  are  given  a  system  whose  state  space  fo  =  A 
X  B,  is  represented  as  the  Cartesian  product  of  two  subsystem 
state  spaces,  A  and  B.  Suppose  further  that  x  ^  si  is 
distributed  according  to  the  pdf  P(x).  If  we  write  x  =  (u,v), 
where  u  ^  A  and  v  B,  then  we  can  obtain  in  the  usual  fashion 
the  marginal  probabilities  of  u  and  v  as 

P.(u)  =  f  P(x)  dv 
A  j 

and  ^  (2) 

P„(v)  =  r  P(x)  du 
B 

We  can  now  obtain  a  second  pdf  on  ft  as  follows: 

Q(x)  =  Pg(v)  . 
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The  random  vectors  u  and  v  are  independent  if  and  only  if 
Pfx)  =  Q(x),  by  definition. 

Now,  given  any  two  pdf’s  Pj^  and  Pg  on  a  single  state 
space,  ft,  J.W.  Gibbs  has  defined  the  function, 

G{Pj,P2)  =  SUM^  C  Pj^j  logtPj^^/Pg^)  ]  ,  (3) 

and  has  proved  that  it  is  always  nonnegative,  and  that  it 
vanishes  if  and  only  if  P^^  =  Pg^  for  all  i  .  It  fails  to  be 
a  metric  on  the  space  of  pdf's  on  ft,  in  part  because  it  is  not 
symmetric  Calthough  that  is  easily  remedied),  but  we  need  not  go 
into  much  more  detail  than  this. 

In  the  special  case  in  which  Pg  is  derived  from  P^^  as  Q 
was  derived  from  P  above,  we  can  now  define  the  structure 
function  Jp(A,B): 

Jp(A,B)  =  G(P.Q)  . 

The  notation  on  the  left  stresses  the  fact  that  the  state  of  the 
original  system  ft  is  distributed  by  the  pdf  P,  which  is  the 
only  pdf  in  sight,  and  that  each  partitioning  of  the  system  into 
factor  spaces  A  and  B  produces  the  nonnegative  number  Jp(A,B) 
which  depends  on  P  and  on  the  two  marginal  probabilities  that 
fall  out  of  the  state  space  factorization.  These  two  marginal 
probabilities  have  their  own  entropies, 
can  further  be  shown  that 

Jp(A,B)  =  E(P^)  +  E(Pg)  -  E(P)  >  0  .  (4) 

The  proof  is  in  Watanabe  [10].  Because  of  this  equation,  the 
structure  function  is  also  referred  to  as  the  "excess  entropy" 
generated  by  the  assembly  of  two  systems  into  one. 
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Note  that  the  structure  function  vanishes  if  and  only  if  the 
subsystems  are  statistically  independent,  i.e.,  if  and  only  if 
their  joint  pdf  (P)  is  the  product  of  their  separate  pdf's.  This 
in  turn  implies  that  there  is  no  (pairwise)  correlation  between 
any  component  of  u  €  A  and  any  component  of  v  €  B,  although 
there  may  well  be  internal  correlations  within  u  and  within  v. 
The  converse  is  false  (lack  of  pairwise  correlation  does  not 
imply  independence),  but  the  contrapositive  is,  of  course,  true: 
If  crosscorrelations  are  nonzero,  then  the  structure  function 
will  be  strictly  positive. 

We  have  defined  a  normalized  version  of  the  structure 
function,  called  DMORPH,  which  is  constrained  to  take  values 
between  0  and  1,  regardless  of  the  dimensions  of  the  state  spaces 
and  regardless  of  the  (finite)  number  of  primitive  events  which 
partition  the  state  spaces.  Thus,  we  have 

DMORPH  =  Jp(A,B)/Mp(A,B),  (5) 


where 


Mp(A,B)  =  E(P)  -  MaxCE(P^),E(Pg)3  . 

The  normalization  divisor,  Mp(A,B),  is  only  our  best  current 
estimate  of  an  upper  bound  for  the  structure  function.  We  have 
not  proven  that  is  a  supremum.  It  is  obtained  by  supposing  that 
the  smaller  of  the  two  systems  (say,  B)  is  completely  correlated 
with  a  subsystem  of  the  larger,  in  which  case  the  "excess 
entropy"  is  the  difference  between  the  entropy  of  the  whole  eund 
the  entropy  of  the  larger  subsystem.  This  un-rigorous  argument 
is  supported  by  numerical  experimentation  and  might  be  made  into 
a  proof  with  a  little  more  effort. 

The  normalized  structure  function,  DMORPH,  is  the  tool  which 
we  have  used  in  our  experiments  to  measure  the  relative  similari¬ 
ty  of  structure  between  two  systems  based  on  vector  samples  of 
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signals  from  the  two  systems.  Computation  of  DMORPH  relies  on 

the  ability  to  estimate  the  three  constituent  entropies,  and  this  j 

is  not  an  easy  task.  Perhaps  the  most  significarit  contribution 
of  this  work  is  the  development  and  implementation  of  this 
algorithm. 

i 

2.3  Computation  of  the  Structural  Similarity  (DMORPH) 

In  order  to  estimate  the  entropy  (relative  or  otherwise)  of 
a  system,  it  is  necessary  to  obtain  an  estimate  of  the  probabil-  ! 

ity  density  function  for  the  state-vector  of  the  system.  In 
order  to  account  for  both  spatial  and  temporal  correlations  in 
the  joint  PDF,  the  state-vector  must  be  sampled  and  held  over  a 

time  interval  which  is  long  enough  to  span  the  coherence  of  the  I 

system.  Since  the  long-term  memory  of  a  neural  network  is 
supposed  to  maintain  temporal  coherence  over  the  lifespan  of  the 
network,  it  will  clearly  not  be  possible  to  sample,  hold,  and 
process  the  full  quantity  of  data  needed  to  characterize  the 
system! 

Instead,  it  will  have  to  suffice  to  use  a  time  window  which 
is  long  enough  to  cover  the  short-term  dynamics  of  the  network, 
i.e.,  its  "impulse  response".  With  such  a  window,  the  resulting 
estimated  PDF  will  reflect  the  short-term  memory  (STM)  of  the 
network,  including  both  the  neuronal  treunsitions  and  the  network 
effects,  but  ceui  only  hope  to  reflect  as  much  of  the  long-term 
memory  dynamics  as  are  evident  within  the  sample  window.  These 
may  not  be  insignificant.  According  to  Klopf  [61  the  coherence 
needed  to  obtain  storage  of  long-term  memories  in  animals  is  on 
the  order  of  about  three  seconds. 

Consider  now  the  problem  of  estimating  the  PDF  for  the 
signal  vector  (including  a  number  of  time  samples)  of  the  system. 

Traditionally,  this  would  be  done  by  partitioning  the  ranges  of 
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the  samples  into  "event  bins",  and  accumulating  a  histogram. 

The  entropy  would  then  be  estimated  by  computing 

E  =  SUM. C(N. /Ni*log(N. /N)],  (6) 

where  ‘L  is  the  number  of  occurrences  of  the  i-th  event,  and  N  is 
the  total  number  of  trials  in  the  experiment.  For  a  vector  with, 
say,  8  components  and  a  time-window  of,  say,  32  samples  per  trial 
(10  samples/sec  for  3.2  sec)  and  a  partition  of  the  range  of  each 
of  those  256  variables  into,  say,  8  levels  —  that  comes  to  n  = 
8**256  possible  events!  This  is  clearly  beyond  the  capacity  of 
available  computational  methods. 

The  estimation  technique  known  as  the  Maximum  Entropy  Method 
(MAXENT,  cf . ,  Jaynes  C5])  overcomes  these  problems  by  constrain¬ 
ing  the  pdf’s  of  interest  to  lie  within  certain  limited  classes 
of  functions.  For  example,  they  may  seek  the  pdf  of  maximum 
entropy  among  all  those  pdf’s  whose  mean  and  covariance  are  equal 
to  the  sample  mean  and  the  sample  covariance.  Under  these  and 
similar  constraints,  the  solution  may  be  found  by  the  method  of 
Lagreurjge  multipliers.  Although  these  methods  have  been  applied 
to  the  study  of  living  neural  networks  [9]  we  assert  that  the 
effort  is  futile.  The  pdf  of  the  signal  vector  of  a  neural 
network  is  typically  extremely  complex,  as  befits  a  system  which 
by  design  extracts  and  stores  millions  of  similarity  clusters  of 
data  which  it  finds  in  its  sensory  inputs.  Thus  the  signals  from 
such  systems  will  of  the  essence  be  "mega-modal".  But  the  MAXENT 
distribution  determined  by  the  first  two  statistical  moments  (on 
an  infinite  domain)  will  be  unimodal,  namely  the  multivariate 
Gaussian.  Since  it  is  precisely  the  fine  structure  of  the  signal 
density  that  we  are  interested  in,  auid  not  the  textbook 
statistical  parameters,  we  need  a  nonparametric  method  which 
reflects  only  those  constraints  imposed  by  our  measuring 
instruments . 
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In  the  end,  our  methods  for  understanding  and  reverse 
engineering  the  brain  will  look  very  much  like  the  method  that 
the  brain  uses  for  underst2unding  and  reverse  engineering  its 
sensory  environment.  Such  an  "eternal  golden  braid"  will  make 
the  recursive  tangles  of  Godel,  Escher,  or  Bach  look  like  mere 
children’s  toys.  In  the  following  description  of  our  method  for 
estimating  structural  similarity,  the  reader  is  invited  to 
observe  that  the  design  of  the  computational  methods  may  hold 
more  glues  to  the  design  of  intelligent  systems  than  will  resu 
from  their  application. 

We  are  approaching  the  problem  in  the  following  way.  For 
computation  of  uhe  DMORPH  function  THREE  entropies  must  be 
computed.  We  can  reduce  that  to  just  ONE  entropy  via  the 
following  argument. 


We  argue  that  ani:  partitioning  of  the  sample  space  into 
events  prior  to  the  collection  of  data  imposes  an  unwarranted 
bias  on  the  resulting  entropy  measurement.  For  example,  the  use 
of  7  threshold  levels  (defining  8  events)  on  a  real-valued 
measurement  must  of  necessity  define  two  regions  which  are 
infinite  in  extent,  and  the  placement  of  these  thresholds 
presumes  some  knowledge  of  where  the  bulk  of  the  measurements 
will  lie.  Therefore,  for  our  computations,  we  use  only  the  a- 
priori  knowledge  of  the  computational  resources  at  our  command  to 
select  the  NUMBER  of  event-bins  for  each  component  of  the  seunpled 
data^  ;  and  we  then  adjust  the  BOUNDARIES  of  these  bins  (i.e., 


1 ;  In  truth,  even  this  strategy  imposes  hidden  a-priori 

constraints.  It  assigns  an  arguably  unwarranted  priority  to 
numerical  contiguity  within  events.  Why,  for  exeimple,  should 
numerically  contiguous  events  be  preferred  by  nature  over, 
say,  a  partitioning  in  which  events  are  defined  by  the  value 
of  the  third  significant  digit  of  an  octal  representation  of 
the  measurements?  The  answer  is  reasonable  and  straight¬ 
forward;  Our  measuring  instruments  have  inertia  and  this 
results  in  an  unavoidable  time-averaging  of  results.  Thus, 
the  definition  of  events  by  contiguous  ranges  of  measurements 
automatically  incorporates  this  smoothing  constraint. 
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the  thresholds  that  separate  them)  during  data  collection  so  that 
the  number  of  observed  events  per  bin  is  the  same  (+/-  1)  for 
each  bin.  This  results  in  a  ‘‘tiling"  of  the  n-dimensional  event 
space  (n  =  product  of  sample-vector  dimension  times  the  number  of 
time-samples  per  trial)  into  equiprobable  events. 


The  actual  procedure  is  illustrated  for  a  two-dimensional 
random  vector  as  follows.  (See  Figure  2. )  Tiling  of  the  subsys¬ 
tem  sample  spaces  is  easily  accomplished  using  a  sort  routine  on 
the  components  of  the  sample  vector.  We  are  currently  partition¬ 
ing  each  component  into  only  two  equiprobable  events  because  of 
the  computational  limitations.  (There  is  a  big  difference 
between  (2)**256  events  and  (3)*»256  events.)  Thus  we  loolt 
first  at  component  i=l  of  the  sampled  data  and  use  a  sort 
routine  to  find  the  median,  where  we  locate  the  single  threshold 
for  the  first  component.  Then,  for  all  sample  vectors  whose 
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first  component  lies  below  the  threshold,  we  find  the  median 
value  of  their  second  components.  This  becomes  the  first  of  TWO 
thresholds  for  the  second  component.  The  second  threshold  is 
obtained  by  finding  the  median  of  the  values  of  the  second 
components  whose  first  components  lay  above  the  threshold.  This 
results  in  four  "tiles"  which  partition  the  two-dimensional 
samples  into  equiprobable  events.  If  the  sample  vectors  were 
three-dimensional,  there  would  be  four  thresholds  on  the  third 
axis,  and  the  seven  total  thresholds  would  partition  the  three- 
dimensional  space  into  8  equiprobable  events. 


Figure  3.  Tiling  of  2-D  Sample  Space  with  8  Events  per  Segment 


This  tiling  of  the  space  is  (in  the  limit  of  large  numbers 
of  thresholds  on  each  axis)  equivalent  to  obtaining  a  PDF  for  the 
data,  in  that  the  reciprocal  of  the  volume  of  each  nonvacuous 
tile  is  proportional  to  the  probability  of  finding  the  state  of 
the  system  in  each  unit  volume  within  that  tile.  Although  the 
tiling  is  obviously  dependent  on  the  order  in  which  the  axes  are 
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selected,  we  conjecture  that  the  differences  become  insignificant 
as  the  number  of  thresholds  on  each  eixis  increases.  It  is 
instructive  to  observe  that  the  tiling  is  a  constructive  repre¬ 
sentation  of  the  probability  measure  which  has  the  histogram  PDF 
for  its  Radon-Nikodym  derivative.  The  measure  (more  precisely, 
the  inner  Jordan  contents  of  any  set  of  states  is  obtained  by 
counting  the  number  of  tiles  within  the  set,  the  seune  as  would  be 
done  if  the  tiles  were  defined  more  customarily  as  unit  hyper¬ 
cubes.  This  is  more  easily  seen  when  the  number  of  events  per 
segment  is  larger  than  2,  as  in  Figure  3.  That  is,  whereas  the 
unit  tiling  which  is  used  for  building  a  normal  histogram  repre¬ 
sents  the  (translation  invariant  Lebesgue)  measure  associated 
with  the  uniform  PDF,  the  tiling  derived  from  the  data  represents 
the  (in  general,  non-translation  invariant)  measure  associated 
with  the  actual  PDF  of  the  data. 

Now,  if  we  had  only  ONE  entropy  to  compute,  there  would  be 
nothing  to  it,  because  our  tiling  guarantees  that  the  sample 
frequencies  are  uniform,  and  the  entropy  of  a  uniform  distribu¬ 
tion  over  n  events  is  the  maximum  possible:  log(n).  But  we  are 
measuring  the  entropies  of  two  presumably  coupled  systems,  having 
n  and  m  events,  respectively,  and  there  are  THREE  entropies  in 
question,  namely  the  two  entropies  of  the  separate  systems,  and 
the  entropy  of  the  composite  system.  We  are  at  liberty  to  tile 
the  sample  spaces  of  the  two  subsystems  any  way  we  like.  But 
once  the  events  in  the  subsystems  are  defined,  then  the  events  in 
the  composite  system  are  determined  as  the  product  space  of  the 
two  subspaces.  The  resulting  entropy  estimate  E(P)  for  the 
composite,  therefore,  must  be  computed  by  the  usual  formula  (6). 
It  will  be  somewhat  less  than  its  meocimum  value  (log(m*n)), 
according  to  the  degree  of  crosscorrelation,  or  structure, 
between  the  subsystems,  and  the  relative  (excess)  entropy  will  be 

Jp  =  log(m)  +  log(n)  -  Elf)  >  0.  (7) 


9 


9 
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From  this,  it  is  then  easy  to  compute  DMORPH  using  equation  (5). 
One  observes  that  because  the  tiling  maximizes  the  entropy 
estimate  for  each  of  the  subsystems,  DMORPH  measures  only  the 
crosscorrelations  which  exist  between  them  and  is  unaffected  by 
any  changes  in  the  internal  organization  of  one  which  are  not 
reflected  by  corresponding  changes  in  the  other.  This  is  also 
confirmed  by  our  experiments. 

Listings  of  the  program  (DTEST)  which  computes  DMORPH  and 
its  associated  tilings  and  entropies  are  included  in  Appendix  B. 
Experiments  showing  the  performance  of  DMORPH  are  described  in 
Section  2.4,  and  the  experimental  conf igurations  and  their 
resulting  graphs  are  shown  in  Appendix  C. 

We  conjecture  that  this  method  of  event-boundary  adjustment 
can  form  the  basis  for  a  learning  law  for  neural  networks  which 
maximizes  the  information  content  of  internal  representations  of 
external  events.  This  will  be  investigated  further  in  Phase  II. 

2.4  DMORPH  Characterization  Experiments 
2.4.1  De.5jCl±Etil2G 

The  DMORPH  experiments  were  performed  to  test  and  verify  the 
performance  of  the  algorithms  which  construct  equiprobable  event 
tilings  of  the  two  subsystem  sample  spaces,  which  compute  the 
entropy  of  the  composite  system,  and  which  compute  the  normalized 
structure  function,  DMORPH. 

I 

In  the  first  set  of  experiments  (1  through  4),  random 

vectors  X  and  Y  of  dimensions  2,  4,  6,  and  8  were  generated. 

2 

Their  components  were  uniform  in  the  interval  [0,1).  These 

I  - 

2:  The  tiling  algorithm  worked  so  well  that  it  immediately 

detected  a  serious  fault  in  our  random  number  generator, 
which  we  subsequently  replaced  with  an  algorithm  from 
Abramowits  &  Stegun. 

I 
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experiments  determined  the  running  time  for  the  algorithm  and 
demonstrated  the  relationship  between  the  dimensions  of  the 
subsystems  and  the  amount  of  data  which  was  needed  to  stabilize 
the  entropies  and  the  structure  value,  DMORPH.  A  sample  plot 
showing  the  entropies  (upper  three  curves)  and  DMORPH  (lower 
curve)  as  they  evolve  with  additional  trials  of  the  experiment  is 
shown  in  Figure  4.  Similar  plots,  together  with  the  experimental 
configurations  are  found  in  Appendix  C  to  document  the 
characterization  experiments. 


EXP36.DAT;  ROUS:  1  TO  200  :  PLOT  OF  TRIAL 


UH.  ENTROPV 
K  ENTROPV 
y  ENIROPy 
-  DMORPH - ^ 


1.0  23.1  45.2 


. . . * . 

- r - • . ■•1--  — ^ 

&7.3  89.4  111.6  133.7  155.8  177.9  200.0 

TRIAL 


Figure  4.  DMORPH  Experiment  for  Dim(X)=2,  Dim(Y)=3 
Pseudorandom  Data  without  gross  correlations. 
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In  the  second  set  of  experiments  (5  and  6),  intracorrela¬ 
tions  were  introduced  into  one  or  both  of  the  random  vectors  X 
and  Y  to  determine  whether  the  measured  structural  similarity 
between  X  and  Y  were  indeed  independent  of  these  intracorre¬ 
lations. 

In  the  third  set  of  experiments  (7  through  35 >,  crosscorre¬ 
lations  were  introduced  between  X  and  Y  to  determine  their 
effect  upon  the  value  of  DMORPH  and  to  evaluate  the  normalization 
divisor  to  see  whether  it  is  close  to  the  theoretical  maximum  of 
the  relative  entropy. 

2.4.2  Results 

The  graphs  in  Appendix  C  illustrate  the  results  of  the 
DMORPH  characterization  experiments. 

Graphs  of  experiments  1  to  4  show  that  as  the  dimensionality 
of  the  systems  increases,  the  quantity  of  data  needed  to  obtain 
stable  estimates  of  the  entropies  and,  hence,  of  DMORPH  also 
increases.  This  is  because  the  total  number  of  distinct  events 

defined  by  the  tiling  algorithm  on  a  subsystem  of  dimension  N 

N 

is  2  ,  and  before  the  subsystems  can  possibly  exhibit  their 
maximum  entropy,  the  number  of  samples  per  bin  must  be  somewhat 
larger  than  unity.  The  tiling  algorithm  insures  that  within  each 
subsystem,  the  events  will  be  defined  so  that  no  matter  how  much 
data  is  taken,  the  greatest  difference  between  the  number  of 
samples  assigned  to  any  two  event-bins  will  be  plus  or  minus  one. 
(Exceptions  occur  whenever  certain  sample  values  occur  multiple 
times,  which  prevents  insertion  of  a  threshold  to  separate  them. ) 
But  whenever  that  difference  is  a  sizeable  fraction  of  the  total 
number  of  data  samples  per  bin,  the  entropy  will  be  substantially 
below  its  theoretical  maximum.  In  particular,  when  the  first 
sample  is  taken,  the  distribution  of  samples  in  the  event  bins 
looks  like  a  delta  function,  so  the  entropy  always  starts  at 
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zero.  After  that,  it  climbs  toward  its  theoretical  maximum 

(since  the  event-bin  boundaries  are  being  adjusted  during  data 

N 

collection),  which  is  log(2  ).  We  convert  the  logarithms  to  the 
base  2,  so  that  the  theoretical  maximum  entropy  for  any  of  these 
systems  is  just  the  dimension  of  the  system. 

Note  that  the  value  of  DMORPH  in  these  graphs  begins  at  zero 
and  then  rises  to  a  maximum  value  before  falling  back  toward 
zero.  The  reason  that  it  falls  back  toward  zero,  of  course,  is 
because  there  is  very  little  crosscorrelation  between  the 
components  of  X  and  of  Y.  (The  residual  makes  a  nice  measure 
of  the  quality  of  the  pseudorandom  number  generator. )  But  the 
fact  that  it  rises  to  a  maximum  showing  some  "false"  structure 
before  adequate  data  is  collected  (see  Figure  4  again)  leads  to 
some  interesting  comparisons  with  the  way  people  learn  about 
their  environment.  It  seems  to  say  that  when  we  are  confronted 
with  a  totally  structureless  system  to  observe,  and  we  begin  to 
collect  data  on  it,  we  will  first  be  convinced  that  the  system  is 
structureless;  then  we  will  begin  to  see  patterns;  but  as  the 
data  becomes  statistically  complete  all  the  patterns  disappear. 

It  also  confirms  the  wisdom  of  carving  out  low-dimensional 
analytical  tasks,  because  with  the  really  big  problems  (e.g.. 
Neural  Network  theory,  or  Artificial  Intelligence)  the  amount  of 
data  which  one  is  likely  to  obtain  during  the  attention-span  of 
the  average  funding  eigency  will  most  likely  lead  one  to  make 
grandiose  claims  of  great  discoveries  which  are  doomed  to 
evaporate  when  the  data  are  more  complete. 

Subsequent  experiments,  described  below,  will  show  the  same 
qualitative  behavior,  except  that  when  there  really  is  some 
crosscorrelation  between  the  subsystems,  the  "false"  structure 
then  gives  way  to  a  nonzero  asymptotic  value. 

In  the  next  set  of  experiments,  numbered  5  and  6  in  Appendix 
C,  intracorrelations  were  introduced  (e.g.,  :=  Xg)  and  the 

results  show  that  DMORPH  is  impervious  to  such  deception. 
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Finally,  crosscorrelations  were  introduced  (e.g.,  Xg  :=  Ygi 
and  Xg  ■■  =  these  show  that  the  structure  function  is 

properly  sensitive  to  them.  Many  experiments  were  performed, 
showing,  e.g.,  the  structure  between  the  input  random  vector  and 
the  output  random  vector  for  various  simple  matrix  transforms. 
Finally,  when  Y  is  made  equal  to  a  subvector  of  X,  DMORPH 
rises  almost  to  unity,  showing  that  the  normalization  divisor  is, 
if  not  precisely  correct,  quite  adequate  to  measure  the  relative 
structure  between  two  systems  without  being  biased  by  the 
dimensionality  of  the  problem. 
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3.  ADAPTIVE  STRUCTURE  OPTIMIZATION 

In  order  to  test  the  applicability  of  the  structure  measure 
to  the  control  of  neural  network  designs,  it  is  necessary  to 
implement  one  or  more  typical  neural  network  designs  and  simulate 
the  interaction  of  that  network  with  the  biological  network. 

With  such  a  simulated  interaction,  it  is  possible  to  treat  the 
randomly  generated  (and  naturally  generated)  input  signals  as 
being  representative  of  the  natural  network  structure,  while  the 
output  signals  represent  the  artificial  network  structure. 

Originally,  it  was  our  intent  to  implement  these  trial 
networks  on  the  MassComp  MC5700  at  the  Biosciences  Laboratory  at 
North  Texas  State  University,  which  has  direct  access  to  the 
signals  emanating  from  a  living  culture  of  several  hundred 
mammalian  neurons  (see  Appendix  A).  However,  it  would  not  be 
possible  to  perform  the  pulse-rate  demodulation  and  the  DMORPH 
tiling  in  real  time,  and  there  is  as  yet  no  capability  for  the 
artificial  network  model  to  talk  back  to  the  natural  network  with 
any  form  of  stimulation.  Therefore,  we  opted  to  take  simulta¬ 
neously  recorded  data  from  multiple  channels  in  digital  form  and 
to  perform  the  experiments  off  line  at  our  own  facilities. 

In  the  following  sub-sections,  we  describe  the  models  which 

we  have  implemented  for  these  experiments.  The  details  of  the 

experiments  2ind  their  results  are  described  subsequently  in 

Section  6  of  the  report.  All  of  our  network  model  software  was 

implemented  under  our  proprietary  dynamical  system  simulator 
TM 

package,  SYSPRO  .  This  has  allowed  us  to  program  the  published 
versions  of  these  models  into  a  flexible  simulation  module  with 
minimal  duplication  of  effort.  In  all  cases,  it  is  only 
necessary  to  produce  a  SYSPRO  primitive  system  which  computes  the 
transfer  function  and  the  learning  algorithm  of  the  subject  model 
and  to  link  it  into  a  possibly  minor  modification  of  our  SYSPRO 
composite  network  model  as  the  replicated  node.  The  specific 
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interconnect  graph  is  specified  at  run  time  through  the  system 
initialization  instructions. 


3.1  The  Back-Propagation  Model 

The  back-propagation  model  which  we  used  is  the  one  which  is 
described  in  Rumelhart  and  McClelland  [7].  This  required  almost 
no  programming  effort,  since  it  is  a  model  with  which  we  have 
extensive  experience  and  which  we  include  as  a  sample  network  in 
our  commercial  neural  network  simulator  package.  The  network  was 
configured  as  a  4-3-4  feedforward  network  (including  direct  links 
from  the  input  layer  to  the  output  layer). 

The  transfer  equations  for  the  processing  elements  are  given 
by 


y  .  =  (r  (  SUM.  z  ..  x.  ;  B.  C  ). 

J  1  0 1  1 

where  «r(A;B,C)  is  a  sigmoidal  function  of  the  first  argument, 
with  values  ranging  between  0  and  1,  whose  maximum  slope,  C,  is 
attained  at  A  =  B.  Except  for  the  ability  to  control  the  value 
of  the  slope  (C)  at  the  desired  threshold  (B),  this  function  is 
the  commonly  used  "logistic"  function: 

<r(x;b,c)  =  l/[  1  +  exp(-4c{x-b)  ]. 

The  learning  law  is  modified  only  so  that  those  factors  of  the 
update  equations  which  can  be  computed  "locally"  are  computed 
'  -thin  the  neuron  model,  and  those  which  require  information  at 
the  network  level  are  computed  by  the  network  model.  (Our 
simulator  protocol  facilitates  the  assembly  of  system  models 
written  in  heirarchical  fashion,  but  it  imposes  the  discipline  of 
using  only  data  which  is  available  to  subsystems  through  the 
input  terminals  and  the  local  state  vector. ) 
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3.2  The  Bear-Cooper-Ebner  Model 

The  BCE  model  is  based  on  the  description  by  Bear,  Cooper, 
and  Ebner  Cl]  of  a  learning  algorithm  attributed  to  Bienenstock. 

This  learning  algorithm  is  partly  Hebbian  and  partly  anti-Hebbian 
in  that  each  synaptic  weight  learns  in  proportion  to  the  presyn- 
aptic  activation,  but  the  proportionality  constant  may  be  posi¬ 
tive  or  negative  according  to  the  relationship  of  the  current 
post-synaptic  activity  to  the  recent  average  of  the  post-synaptic 
activity.  That  is.  if  the  recent  activity  has  been  high,  but  the 
current  activity  is  lower  than  the  average  activity,  the  synaptic 
weight  will  be  reduced.  If  the  recent  activity  has  been  low, 
then  almost  any  post-synaptic  activity  will  be  greater  than  the 
average  activity  and  will  cause  an  increase  in  the  synaptic 
weight.  The  Bienenstock  law  is, 

dm  ./dt  =  0(c, c)  d  .  , 

0  J 

( 

where  m.  is  the  j-th  synaptic  weight,  d.  is  its  presynaptic 
signal,  c  is  the  neuronal  output  signal  (in  the  linear  region), 
and  c  is  the  average  of  c  over  a  recent  time  interval. 

We  implemented  the  0  function  (see  Figure  5)  as  a  spline 
of  a  parabola  on  the  left  and  an  exponential  learning  curve  on 
the  right  of  the  crossover  point,  0j^  .  We  implemented  the 
neuronal  transfer  function  more  generally  than  is  described  in 
BCE,  so  as  to  include  a  sigmoid  nonlinearity  at  the  output  (the 
same  logistic  sigmoid  used  above  in  the  back-propagation  model), 
rather  than  to  assume  operation  in  the  linear  region.  This 
necessitated  a  decision  on  the  interpretation  of  c  in  the 
learning  law,  and  we  chose  to  interpret  c  as  the  neuronal 
output,  rather  than  as  the  postsynaptic  activation. 
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Figure  5.  The  Bienenstock  0  Function  (two  versions). 


3.3  The  Klopf  "Drive-Reinforcement"  Model 

We  implemented  the  drive-reinforcement  model  in  accordance 

with  the  description  given  in  Klopf  [63.  In  order  to  control  the 

learning  rate,  we  augmented  the  learning  algorithm  with  a  GAIN 

factor,  which  in  effect  scales  the  area  under  the  learning  rate 

constant  curve  (the  curve  determined  by  the  constants,  c  .,  in 

Klopf’ s  report).  By  setting  GAIN  =  0,  we  could  turn  learning  off 

at  any  time  so  that  the  tiling  operation  of  DMORPH  would  have  a 

time- invariant  segment  of  the  network’s  signals  to  work  with, 

i.e.,  one  in  which  the  structure  was  not  changing  during  the 

tiling  operation.  Also,  since  the  GAIN  factor  was  included,  we 

chose  to  implement  the  learning  rate  constants,  c^,  so  that 

their  sum  is  unity.  This  then  treats  the  c  .  vector  as  a  weight 

0 

vector  for  a  weighted  average  of  the  prior  history  of  the 
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synaptic  efficacies.  The  relative  magnitudes  of  our  values  for 
c^  were  (almost)  the  same  as  used  by  Klopf. 

When  the  Klopf  neurons  were  connected  into  a  network,  we 
decided  not  to  make  any  of  the  synapses  non-plastic,  since  inside 
the  network  it  is  unlikely  that  the  distinction  between  condi¬ 
tioned  stimuli  (CS)  and  unconditioned  stimuli  (US)  could  be  made 
a-priori,  and  in  any  case  our  time  constraints  did  not  allow  such 
fine-tuning. 
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-  This  pa^e  is  mostly  blank 
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4.  DESIGN  OF  THE  BIDIRECTIONAL  AXON  BUNDLE 

In  this  section  we  describe  the  proposed  technique  for 
constructing  an  interface  between  a  living  tissue  culture  of 
active  mammalian  neurons  and  an  artificial  neural  network  which 
is  hosted  in  a  general  purpose  computer.  The  design  is  based  on 
the  laboratory  setup  in  the  neurophysiology  laboratory  of  Dr. 
Guenter  W.  Gross  at  North  Texas  State  University  and  on  his 
proposed  apparatus  for  localized  stimulation  of  that  network. 

First  we  present  a  brief  description  of  the  NTSU  laboratory 
apparatus,  which  is  more  thoroughly  described  in  the  Appendix. 
After  that,  we  describe  the  status  of  the  work  being  conducted  at 
NTSU  and  at  Southern  Methodist  University  to  achieve  the  local¬ 
ized  stimulation  of  the  culture  network.  Finally,  in  the  third 
subsection  following,  we  describe  the  functional  design  of  a 
bidirectional  interface,  called  the  Synthetic  Axon  Bundle,  which 
will  enable  the  culture  network  to  influence  and  be  influenced  by 
the  signals  in  an  artificial  neural  network. 


4. 1  Description  of  the  NTSU-Gross  Apparatus 

Professor  Gross’s  laboratory  apparatus  is  described  briefly 
here  and  illustrated  thoroughly  in  Appendix  A.  The  multimicro¬ 
electrode  plate  (MMEP)  on  which  the  culture  is  maintained  is 
described  first,  followed  by  a  description  of  the  recording 
chamber  design.  The  digital  processing  system  is  illustrated  on 
page  A-18. 

Signals  from  the  neural  culture  are  amplified  and  patched 
into  the  data  acquisition  and  control  processor  of  the  Masscomp 
MC5700  computer,  where  they  are  converted  to  digital  form, 
typically  at  a  30  kHz  per  channel  sample  rate.  The  digital 
signal  is  filtered  for  A/C  hum  and  is  then  processed  for  burst 
detection  (pages  A-24  to  A-29).  Signals  at  various  stages  of 
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processing  can  be  selected  for  display  on  the  color  monitor  (page 
A-25). 


The  signal  monitoring  apparatus  is  supplemented  by  an 
auditory  monitor,  which  codes  each  electrode's  activity  into  tone 
bursts  at  a  frequency  that  is  unique  to  the  source  electrode,  and 
by  an  LED  display  which  provides  visual  cues  to  the  activity  on 
each  electrode.  These  were  originally  designed  as  "PR  enhance¬ 
ment  instruments"  (where  PR  stands  for  Public  Relations),  but 
they  have  proven  to  be  valuable  intuitive  aids.  The  human  ear 
can  detect  patterns  and  correlations  in  the  data  that  would  go 
completely  unnoticed  on  a  strip-chart  recording. 

Not  shown  in  the  hardware  description  of  Appendix  A  is  a 
limited  capability  for  stimulation  of  the  cultured  network.  This 
is  described  in  greater  detail  below.  First,  we  describe  the 
design  of  a  Synthetic  axon  bundle  which  presumes  the  availability 
of  a  suitable  stimulation  apparatus. 


4.2  Functional  Design  of  the  Synthetic  Axon  Bundle 

The  purpose  of  the  Synthetic  Axon  Bundle  is  to  provide  the 
communication  link  between  the  natural  mouse  neural  network  (MNN) 
in  culture  and  the  Artificial  Neural  System  (ANS>  oeing  simulated 
on  the  MASSCOMP.  That  is,  it's  function  is  to 

(1)  Modulate  the  signals  that  are  output  from  the  ANS  so 
that  they  can  be  used  to  stimulate  the  MNN  and 

(2)  Demodulate  the  signals  that  are  recorded  from  the  MNN  so 
that  they  can  be  used  as  input  to  the  ANS. 

Figure  6  illustrates  the  basic  functional  design,  and  the 
following  description  provides  some  of  the  details. 
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1 


Figure  6.  Design  of  the  Synthetic  Axon  Bundle 


ANS - >  MNN  ^ 

It  is  believed  that  the  output  of  an  ANS  should  represent 
some  sort  of  information  transfer,  either  by  its  effect  on  the 
modification  of  synaptic  weights  or  the  interpretation  of  what  ft 

this  neuron  firing  means  (the  recognition  of  some  pattern,  say). 

Therefore  it  is  necessary  to  modulate  each  ANS  output  into  a 

spike  train  to  be  used  to  stimulate  the  neurons  of  the  MNN  which 

are  in  close  proximity  to  a  particular  electrode.  There  is  no  ^ 

requirement  that  the  processing  clocks  of  the  two  neural  systems 

be  on  the  same  time  scale.  The  only  requirement  that  exists  is 

that  the  data  from  the  ANS  be  modulated  in  such  a  manner  that  the 

MNN  finds  it  to  be  "stimulating".  Some  experimentation  will  be 
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needed  early  in  Phase  II  to  verify  that  the  following  proposed 
design  will  result  in  a  real  time  stimulus  of  the  MNN. 

The  simulation  time  increment,  ^t,  of  the  ANS  model  will  be 
set  so  that  the  data  produced  by  it  each  At  can  be  modulated  and 
used  as  a  stimulus  over  the  next  At  seconds  for  the  MNN.  Of 
necessity,  there  will  be  a  At  second  time  lag  between  the  data 
output  from  the  ANS  and  the  stimulus  being  applied  to  the  elec¬ 
trodes  of  the  MNN.  Since  we  expect  to  perform  temporal  sampling 
as  well  as  sampling  across  the  electrodes  this  should  not  be  a 
problem  in  our  search  for  cross-system  structure. 

The  amplitude  of  the  voltage  spikes  generated  to  drive  the 
MNN  will  be  consistent  with  the  amplitudes  observed  in  the  MNN. 
The  frequency  used  will  be  proportional  to  the  signal  amplitude 
out  of  the  ANS. 


MNN - >  ANS 

The  signals  recorded  from  the  MNN  were  collected  on  one  of 
multiple  micro  electrodes  in  Dr.  Gross's  laboratory.  The  data 
used  in  many  of  the  experiments  performed  during  the  Phase  I 
effort  was  compressed  using  the  following  processing  algorithm. 


Eight  simultaneous  channels  of  data  were  collected  at  a 
data  rate  of  30,000  Hz. 

This  30,000  Hz.  data  rate  is  then  reduced  to  500  Hz.  by 
saving  only  the  maximum  absolute  value  of  each  disjoint 
and  contiguous  set  of  60  data  values  on  each  channel. 

The  dynamic  range  of  these  data  values  is  further  reduced 
by  comparing  the  data  value  with  a  threshold  and  replacing 
it  by  a  1  if  the  data  value  is  greater  than  or  equal  to 
the  threshold  or  replacing  the  data  value  by  zero  if  it  is 
less  than  the  threshold. 

Finally,  a  16  point  rectangular  filter  (with  unit  weights) 
is  applied  to  a  sliding  window  of  the  data  so  that  the 
data  which  is  input  to  the  ANS  is  an  integer  between  0  and 
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16,  inclusive.  This  data  is  to  be  interpreted  as  Pulse 
Repetition  Frequencies  (PRF)  for  the  neurons  which  are 
being  recorded  on  each  electrode.  The  following  table 
gives  the  range  of  PRFs  for  each  of  the  possible  17  data 
values . 


M(t)  =0 

PRF 

< 

31.25 

Hz 

M(t)  =1 

31.25 

PRF 

< 

62.50 

Hz 

M(t)  =2 

62.5 

PRF 

< 

93.75 

Hz 

M(t)  =3 

93.  75 

PRF 

< 

125.  0 

Hz 

M(t)  =4 

125. 50 

PRF 

< 

156.25 

Hz 

M{t)  =5 

156.25 

% 

PRF 

< 

187.50 

Hz 

M(t)  =6 

187. 50 

,«• 

PRF 

< 

218. 75 

Hz 

M(t)  =7 

218.75 

PRF 

< 

250. 00 

Hz 

M(t)  =8 

250.  0 

PRF 

< 

281.25 

Hz 

M(t)  =9 

281.25 

PRF 

< 

312. 50 

Hz 

M(t)  =10 

312.50 

I?. 

PRF 

< 

343.75 

Hz 

M(t)  =11 

343. 75 

PRF 

< 

375.00 

Hz 

M(t;  =12 

375. 00 

PRF 

< 

406. 25 

Hz 

M(t)  =13 

406.25 

PRF 

< 

437.50 

Hz 

M(t)  =14 

437.50 

PRF 

< 

468.75 

Hz 

M(t)  =15 

468. 75 

PRF 

< 

500. 00 

Hz 

M{t)  =16 

PRF 

> 

500. 00 

Hz 

This  algorithm  has  several  defects,  but  it  also  has  the 
important  advantage  of  its  mere  existence.  Thus,  we  were  at 
least  able  to  run  experiments  on  genuine  digitized  data  from  the 
MMEP,  but  the  interpretation  of  results  must  be  qualified  by  the 
effect  of  the  following  problems. 


First,  the  data  collected  on  any  single  micro  electrode  is 
known  to  be  the  result  of  firings  of  multiple  neurons.  These 
signals  should  really  be  separated  rather  than  lumped  together. 


Second,  the  threshhold  used  to  declare  a  signal  present  on 
the  channel  (electrode)  is  not  changing  over  time  and  therefore 
the  false  alarm  rate  is  not  constant  on  the  channel. 


Third,  the  data  is  collected  at  an  extremely  high  data  rate 
but  the  high  data  rate  features  of  the  data  are  not  exploited  in 
the  signal  processing  algorithm  at  all.  Why  waste  all  the 
magnetic  storage? 
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Fourth,  the  definition  of  bursting  for  the  channel  is 
limited  to  17  discrete  values  rather  than  taking  on  any  value  on 
the  positive  interval  from  zero  to  the  sampling  rate. 

A  set  of  signal  processing  algorithms  which  addresses  the 
majority  of  these  problems  have  been  formulated  by  Martingale 
Research  Corporation  and  supplied  to  Dr.  Gross,  but  they  have  not 
yet  been  implemented  due  to  lack  of  funds  at  NTSU  for  support  of 
student  programmers.  The  design  specification  for  these  algo¬ 
rithms  were  presented  in  (Dawes  and  Collard  [3])  and  are  repro¬ 
duced  in  Appendix  D. 
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5.  PROGRESS  WITH  STIMULATION  CONDITIONING 

In  order  to  utilize  our  structure  function  for  improvement 
of  neural  network  designs  in  a  real  time  interactive  experiment, 
it  is  required  that  the  natural  and  the  artificial  neural  net¬ 
works  be  connected  for  bidirectional  communications.  Part  of  the 
technology  to  accomplish  that  is  described  in  the  preceding  sec¬ 
tion  and  is  called  the  "Synthetic  Axon  Bundle".  It  is  basically 
little  more  than  a  modulator/demodulator  (MODEM)  which  translates 
the  signals  from  a  form  suitable  to  their  source  to  a  form 
suitable  to  their  destination.  But  the  specific  component  which 
injects  the  signal  into  the  MMEP  has  not  been  specified  or  tested 
yet. 


What  is  needed  are  the  following  capabilities: 

1.  A  multichannel  pulse  generator  whose  signal  output 
characteristics  are  subject  to  computerized  control 
individually  by  channel  according  to  pulse  eimplitude  and 
pulse  rate. 

2.  A  localization  of  the  voltage  gradients  within  the  MMEP 
so  that  gradients  capable  of  inducing  depolarization 
into  neurites  are  limited  to  the  vicinity  of  the  active 
electrode. 

The  first  requirement  is  not  a  big  problem,  but  the  second  is 
more  difficult  to  satisfy.  At  present,  the  input  signal  is 
applied  between  the  selected  electrode  and  the  metallic  bezel 
which  surrounds  the  recording  area  (see  page  A-9).  This  results 
in  depolarization  of  an  estimated  10%  to  40%  of  the  neurons  in 
the  culture.  Anything  which  is  more  selective  will  have  to  be 
based  on  the  bipolar  excitement  of  adjacent  pairs  of  MMEP  elec¬ 
trodes,  and  this  will  require  some  redesign  of  the  preamplifier 
boards . 


« 


« 
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A  team  of  electrical  engineers  under  the  direction  of  Prof. 
Lorn  Howard  at  the  Electrical  Engineering  Department  of  Southern 
Methodist  University  is  working  on  the  stimulation  problem.  Dr. 
Gross  at  NTSU  presently  has  a  spike-signal  generator  connected  to 
the  MMEP  which  is  capable  of  injecting  a  signal  into  a  single 
electrode  at  a  selectable  pulse  width  and  rate.  But  without 
sufficient  control  to  be  able  to  simulate  bursting,  he  is  unable 
to  demonstrate  any  form  of  conditioning  of  the  network. 
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6.  EXPERIMENTS  PERFORMED  AND  THEIR  RESULTS 

6.1  Back-Propagation  Experiments 

We  ran  some  simple  experiments  using  an  eleven  neuron  feed¬ 
forward  network  which  learns  by  the  "back-propagation"  algorithm 
to  simulate  the  communication  of  a  noncontrol lable  neural  network 
with  a  controllable  one.  The  feedforward  network  is  the  4-3-4 
network  shown  in  Figure  7.  Its  four  input  signals  consisted  of 
raised  sinusoids  of  various  amplitudes  and  phases  and  its  output 
is  determined  by  the  weights  and  biases  of  the  processing 
elements . 


To  test  the  structure  function  on  this  network,  we  ran  two 
experiments.  Both  experiments  used  the  same  input  vector 
consisting  of  the  four  raised  sinusoids,  and  both  began  with  a 
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random  initialization  of  the  weights  of  the  network.  With 
learning  turned  off,  the  network  was  exercised  for  1500  seconds 
of  simulation  time  (at  10  state  updates  per  second)  to  obtain 
baseline  structure  with  the  randomly  initialized  weights.  Then 
the  learning  was  turned  on  and  the  "desired"  output  was  set  equal 
to  the  input  vector.  Convergence  was  fairly  rapid,  and  the 
output  took  the  qualitative  appearance  of  the  input,  with 
distortions  appropriate  to  the  nonlinearities  of  the  sigmoid 
function.  Then  learning  was  turned  off  and  the  network  was  again 
exercised  to  obtain  the  post-learning  structure. 

In  the  first  experiment  (BP4  1*1,2,  p.  C.  34-35),  the  baseline 
value  of  DMORPH  was  0.356  prior  to  learning  and  0.336  afterward! 
The  structure  actually  declined  after  learning.  Obviously  there 
was  an  error  in  our  procedure. 

Ordinarily,  we  would  not  report  on  the  many  experiments  in 
which  we  have  identified  procedural  errors,  but  in  this  case,  zhe 
error  was  especially  instructive  and  it  supplied  us  with  an 
important  clue  to  the  role  of  causality  in  neural  transductions. 
That  clue  is  taken  up  again  and  discussed  at  greater  length  in 
Section  7.  For  now,  we  merely  point  out  that  in  the  first 
experiment,  we  sampled  both  the  input  vector  and  the  output 
vector  immediately  upon  presentation  of  each  new  input  to  the 
network.  Consequently,  the  signal  that  was  present  on  the  output 
terminals  was  the  one  that  was  left  over  from  the  previous  input. 

The  second  experiment  was  the  same  as  the  first,  except  that 
the  inputs  and  outputs  were  sampled  on  the  half-second  instead  of 
at  each  whole  second  of  the  simulation  clock.  This  allowed  five 
simulation  cycles  for  the  input,  which  only  changes  at  the  begin¬ 
ning  of  each  whole  second,  to  propagate  through  to  the  output. 

In  this  experiment  the  pre-learning  structure  was  0.446  (BP4  <*3, 
Page  C.36),  and  the  post-learning  structure  was  0.515  (BP4  **4, 
Page  C.37). 
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Not  only  did  the  post-learning  structure  show  an  increase 
over  the  pre-learning  structure,  but  the  pre-learning  structure 
as  determined  by  the  lagged  sampling  (0.446>  was  higher  than  the 
pre-learning  structure  of  the  first  experiment  (0.356).  This  is 
because  in  the  first  experiment,  the  structure  reflected  the 
correlation  that  exists  between  one  sample  of  a  sinusoid  and  the 
network-transduced  image  of  another  sample  approximately  1/20 
cycle  into  the  past.  In  the  second  experiment,  the  structure 
reflected  the  correlation  that  exists  between  each  input  sample 
and  its  own  transduced  image. 


6.2  BCE  Experiments 

The  BCE  experiments  were  originally  planned  to  parallel  the 
foregoing  BPE  experiments  for  a  small  network  whose  learning  law 
is  due  to  Bienenstock  (described  above).  After  implementing  and 
testing  a  network  of  neurons  using  that  learning  law,  however,  we 
observed  a  phenomenon  that  led  us  to  discard  our  planned  experi¬ 
ments,  at  least  for  now. 

What  we  observed  was  that  if  a  constant  nonzero  input  was 
supplied  to  one  synapse  of  a  BCE  neuron,  the  synaptic  weight,  and 
consequently  the  neuronal  output,  quickly  reached  a  periodic 
state.  The  period  of  the  oscillation  was  directly  proportional 
to  the  length  of  the  window  used  for  the  sliding-window  average 
of  the  postsynaptic  activation.  That  is,  a  long  window  produced 
a  low  frequency  oscillation,  while  a  short  window  produced  a  high 
frequency  oscillation.  The  length  of  the  period  is  of  the  same 
order  of  magnitude  as  the  length  of  the  sliding  window.  The 
graph  in  Figure  8  shows  this  behavior. 

In  retrospect,  it  might  have  been  worthwhile  to  go  ahead  and 
perform  the  DMORPH  experiments  on  this  kind  of  network,  but  there 
was  not  sufficient  time.  We  report  these  results  as  a  matter  of 
interest  regarding  possible  future  use  of  the  model. 
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Figure  8.  Oscillatory  BCE  Response  to  Unit  Step  Input 


6.3  Drive-Reinforcement  Experiments 

In  late  December,  we  received  Klopf’s  report  [63  on  his  work 
with  the  modified  differential  Hebbian  learning  law  of  Sutton- 
Barto,  which  he  calls  the  Drive-Reinforcement  model.  We  were 
able  to  implement  a  small  network  of  neurons  with  the  D-R 
learning  law  in  a  couple  of  days,  and  we  ran  a  set  of  experiments 
similar  to  the  BPE  experiments,  but  on  a  four-neuron  network.  In 
this  series  of  experiments,  each  neuron  received  a  separate 
input,  and  each  neuron’s  output  was  sampled.  Thus  our  structure 
function  was  evaluated  on  a  four-by-four  system. 
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Figure  9.  The  Drive-Reinforcement  Network 


The  experiments  were  designed  to  measure  the  effects  of  both 
the  learning  algorithm  and  the  architectural  parameters  on  the 
structural  similarity  between  the  input  vector  to  the  network  and 
its  output  vector.  To  determine  the  effect  of  the  learning 
algorithm,  we  began  with  a  random  initialization  of  the  synaptic 
weights  on  a  four- neuron  network  connected  as  in  Figure  9.  This 
interconnection  incorporates  two  feedback  loops,  which  helps  to 
randomize  the  latency  between  onset  of  signals  at  the  various 
inputs  to  each  neuron  and  thus  guarantees  that  the  drive- 
reinforcement  learning  algorithm  sees  the  necessary  delays.  We 
did  not  hard-wire  any  synapses,  since  in  the  network  setting  it 
is  not  clear  that  any  synapse  may  know  a-priori  that  it  will  be 
receiving  the  unconditioned  stimulus  (US).  The  components  of  the 
C  vector  (learning  rate  at  delays  of  0.5,  1.0,  1.5,  2.0,  2.5  sec) 
were  established  at  approximately  1/10  of  the  values  suggested  by 
Klopf  so  that  the  area  under  the  sliding  "C  "  window  would  be 
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unity.  In  particular,  we  chose  C  =  (.4,  .25,  .2,  .1,  .05).  By 
setting  the  GAIN  parameter  at  10,  we  could  scale  these  back  up  to 
approximate  Klopf's  settings,  or  we  could  alter  their  values 
directly. 

The  four  input  signals  were  the  same  raised  sinusoids  as 
were  used  for  the  Back  Propagation  experiments,  except  that  their 
values  only  changed  every  three  seconds  in  this  experiment. 

Since  the  Drive-Reinforcement  neuron  cycles  only  twice  each 
simulation  second,  this  gives  time  for  the  signal  to  propagate 
around  the  feedback  loops  before  the  input  is  changed.  Other¬ 
wise,  the  procedure  was  the  same  as  before.  Learning  was  turned 
off  to  obtain  a  baseline  input/output  structure.  Then  learning 
was  turned  on  until  the  synaptic  weights  had  undergone  signifi¬ 
cant  change  (but  not  long  enough  for  any  weight  saturations  to 
occur).  Then  learning  was  turned  off  and  the  network  was 
exercised  with  the  learned  weights  to  obtain  the  post- learning 
input/output  structure. 

The  pre-learning  structure  with  random  weights  was  0. 480  and 
after  learning  it  was  0.412  (Page  C.38,39). 

A  second  experiment  was  run  in  which  the  values  of  C(l)  and 
C(2)  were  reversed.  This  was  an  attempt  to  determine  the  effect 
of  architectural  adjustments  on  the  post-learning  structure.  The 
pre-learning  structure,  of  course,  did  not  change,  because  this 
architectural  change  does  not  affect  propagation  when  the  learn¬ 
ing  is  turned  off.  After  learning  for  the  same  length  of  time  as 
before,  the  post-learning  structure  was  0.391,  slightly  worse 
than  before  (Paige  C.40). 


6.4  Other  Experiments. 

One  of  our  experiments  was  designed  to  look  into  the  impli¬ 
cations  of  time  and  causality.  That  experiment  was  performed  by 
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partitioning  five  of  the  8  channels  of  actual  action  potentials 
from  the  MMEP  into  an  X  signal  of  dimension  4  and  a  Y  signal 
of  dimension  1.  We  then  extended  the  Y  signal  along  the  time 
axis  so  that  its  first  component  was  sampled  at  the  same  time  as 
all  four  of  X's  components  were  sampled,  and  its  next  three 
components  were  taken  at  three  subsequent  time  values,  so  that  Y 
is  also  a  four-vector.  Our  reasoning  is  that  if  one  or  more  of 
the  X  signals  were  to  correlate  with  anything  in  Y,  it  would 
only  be  seen  at  a  later  time.  The  situation  is  shown  in  Figure 
10. 


"X"  half 
of 

culture 


•’Y‘'  half 
of 

culture 


Y 


Figure  10.  Sampling  Diversity  in  Space  and  in  Time 


The  result  was  a  DMORPH  value  of  0. 1,  which  is  far  less  than  was 
achieved  with  the  simultaneous  sampling  in  which  both  X  and  Y 
were  4-vectors  (DMORPH  =  0.3).  Similar  experiments  on  different 
segments  of  the  MMEP  data  show  qualitatively  similar  results. 

We  also  programmed  a  network  of  Grossberg  Avalanche  neurons, 
which  learn  by  a  form  of  the  basic  Hebb  rule.  We  used  them  to 
verify  the  basic  performance  of  one  of  the  avalanche  architec¬ 
tures,  but  due  to  time  constraints,  we  did  not  employ  these 
models  in  any  of  the  structure  function  experiments.  This  kind 
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of  experiment  on  these  models  will  have  to  await  Phase  II 
research. 

The  results  of  the  initial  tests  have  been  to  illustrate 
that  the  measurement  of  the  relative  entropy  between  two  systems 
is  not  a  simple  matter,  as  one  might  already  guess  from  reading 
the  MAXENT  literature  (of.  [5]).  Nevertheless,  we  have  obtained 
an  algorithm  which  operates  well  within  the  time  and  memory 
constraints  imposed  by  a  PC  computing  environment  when  the  vector 
dimension  of  the  combined  random  signals  from  two  systems  is  8  or 
less.  This  is  adequate  for  determining  the  viability  of  the 
proposed  method. 
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7.  FINDINGS  AND  ANALYSIS 

In  the  Back  Propagation  experiments,  the  structure  increased 
after  learning,  but  in  the  Drive-Reinforcement  experiments,  it 
actually  went  down  afterward.  Moreover,  in  the  MMEP  signal 
studies  we  saw  more  structure  in  the  simultaneous  samples  than  in 
the  samples  showing  both  spatial  and  temporal  diversity. 

These  results  seem  ambiguous  at  best.  In  the  case  of  the 
Back  Propagation  experiments,  the  structure  could  hardly  fail  to 
increase  after  learning,  since  the  "desired  "  output  was  known  to 
be  strongly  correlated  with  the  input  and  the  BPE  network  learns 
to  produce  the  desired  output  rather  well.  Tn  the  Drive- 
Reinforcement  experiments,  the  failure  of  the  structure  to  rise 
after  learning  is  puzzling.  It  is  clearly  an  indicator  that  the 
D-R  learning  law  does  not  necessarily  enhance  correlations 
between  input  and  output  signals,  but  then  it  was  not  designed  to 
do  that,  at  least  not  directly. 

We  devoted  considerable  effort  to  our  attempts  to  understand 
these  results  and  to  appreciate  their  implications  for  both  the 
design  of  learning  laws  and  the  optimization  of  architectural 
parameters.  Eventually,  thanks  to  insistent  questioning  by  our 
student  assistant,  Mr.  David  Boney,  we  happened  upon  the  follow¬ 
ing  "gedankenexperiment" ,  which  shows  very  clearly  that  one 
cannot  naively  infer  that  the  best  neural  network  is  the  one 
which  generates  the  greatest  value  of  DMORPH  between  its  input 
vector  and  its  output  vector. 

Suppose  that  the  artificial  network  were  constructed  with  an 
array  of  "bypass  valves"  corresponding  to  each  of  its  inputs,  and 
that  these  valves  served  to  proportionately  disconnect  the 
network  from  its  inputs  and  reroute  those  inputs  directly  to  the 
output  terminals.  Then,  as  we  turned  the  dials  we  would  see  the 
structural  similarity  between  the  input  signal  vector  and  the 
output  signal  vector  improve  dramatically  toward  its  theoretical 
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maximum,  and  our  inference  would  lead  us  to  conclude  that  the 
best  network  architecture  was  the  one  that  was  not  there  at  all! 

The  problem  is  not  a  fault  with  the  structure  function,  but 
rather  with  its  application  and  the  inferences  drawn  from  it.  We 
hasten  to  point  out  that  what  we  have  done  in  this  project  is  to 
quantify  and  apply  concepts  that  many  neural  network  and  cogni¬ 
tive  science  researchers  have  tacitly  and  Qualitatively  assumed 
to  be  at  work  in  self-organizing  systems.  Our  experiments  have 
shown  that  these  assumptions  need  to  be  much  more  carefully 
thought  out. 

What,  then,  is  the  measure  of  a  cognitive  system?  At  the 
beginning  of  this  work,  we  dismissed  naive  self-organization  — 
at  least  as  it  might  be  measured  by  information-theoretic 
entropy,  since  that  clearly  favors  mental  crystallization.  We 
now  seem  driven  to  dismiss,  or  at  least  to  severely  qualify,  the 
placement  of  any  value  on  the  network's  introduction  of  cross 
correlations  between  its  inputs  and  its  outputs.  Parrots  are 
only  amusing  for  a  short  while,  and  networks  which  merely 
associate  "desired"  outputs  with  selected  clusters  of  inputs 
hardly  know  what  constitutes  a  surprise,  much  less  do  they  have 
any  hope  of  developing  an  appropriate  response  to  one. 

The  defect  in  these  models,  insofar  as  they  attempt  to 
represent  basic  elements  of  cognitive  systems,  is  that  they  pay 
too  little  attention  to  the  fundamental  role  of  time  and 
Gaua.aLi±y.  Even  the  conditioning  models  discussed  by  Klopf  [6], 
which  explicitly  account  for  the  temporal  ordering  of  '^svents  and 
the  temporal  gradients  within  excitations,  and  which  do  an 
admirable  job  of  modeling  drives,  reflexes,  and  even  the  act  of 
generalization,  will  probably  not  emerge  into  cognitive  systems 
through  the  blessings  of  mere  complexity. 
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Figure  11.  The  Neural  System  Design  Problem 


Our  experiments  with  time  and  causality  discussed  above  have 
led  us  to  the  following  description  of  a  learning  system,  which 
pays  special  attention  to  its  relationship  to  its  environment. 

Figure  11  illustrates  the  situation  in  a  manner  which  is  intended 

to  be  especially  significant  to  mathematicians  who  may  have  ^ 

studied  category  theory.  Category  Theorists  are  sometimes  known 

among  mathematicians  as  arrow  chasers".  In  this  case,  we  are 

interested  in  the  fact  that  there  are  three  paths  leading  from 

the  initial  state.  A,  of  the  observed  system  to  the  final  state,  I 

B'/B",  of  the  neural  system. 

One  path,  which  we  call  the  LL  path  (for  Lower  Left), 
represents  the  sensors  mapping  the  initial  state  A  of  the  I 

observed  system  into  an  initial  internal  representation  A’ . 

Then  the  neural  network  transforms  that  representation  into  a 
final  state  B"  (which  may  be  only  an  infinitesimal  time,  dt,  away 
if  we  think  of  these  as  differential  systems). 

I 
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Another  path,  the  UR  path,  represents  the  observed  system 
evolving  from  its  initial  state  A  to  its  final  state  B, 
possibly  under  the  influence  of  inputs  from  the  outside  world. 
Then  the  sensory  network  maps  that  final  state  into  the  internal 
representation,  B' . 

The  third  path,  called  SN  (for  SiNuous),  represents 
observation  of  the  initial  state,  followed  by  action  of  the 
neural  system  upon  the  observed  system  through  motor  controls, 
which  produces  a  controlled  final  state  B.  This  final  state  is 
then  observed,  producing  the  internal  representation,  B' . 

For  the  moment,  we  shall  ignore  the  SN  path,  emd  ask  how  the 
network  can  make  TFDC  (an  acronym  known  by  category  theorists  to 
mean  "The  Following/Foregoing  Diagram  Commutes").  That  is,  how 
can  the  two  representations,  B'  and  B",  be  made  to  coincide,  so 
that  both  the  LL  and  the  UR  paths  produce  the  same  result?  The 
answer  is  that  it  must  build  within  itself  a  state  transition 
operator  which,  when  composed  with  the  effect  of  the  sensors, 
produces  the  same  result.  This  is  the  meaning  of  "learning"  in  a 
sense  which  makes  essential  use  of  the  dynamics  of  the  universe, 
including  that  small  part  of  it  called  the  neural  system.  It 
absolutely  must  employ  a  means  of  comparison  between  the  two 
representations,  B'  and  B".  We  have  illustrated  that  comparison 
by  the  juxtaposition  of  two  state  boxes  for  B'  and  B”,  and  an 
arrow  which  returns  a  control  signal  to  the  current  transition 
operator  of  the  network;  but  we  do  not  mean  to  imply  literally 
that  there  must  be  two  such  "slabs"  within  the  network  together 
with  an  explicit  "metric"  between  them.  That  may  be  the  case, 
but  it  may  also  happen  that  the  sensory  map  and  the  cognitive 
prediction  converge  on  the  same  leo^er  to  produce  a  disturbance 
away  from  homeostatic  equilibrium  at  precisely  those  locations 
harboring  the  pieces  of  the  distributed  transition  operator  which 
need  correcting. 
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It  happens  that  we  have,  in  the  process  of  analyzing  these 
experiments,  constructed  both  a  neural  model  (a  totally  novel  one 
at  that)  AND  an  appropriate  architecture  to  imbed  it  into,  which 
may  achieve  these  goals.  The  model  is  highly  preliminary,  and 
since  it  was  conceived  in  the  final  days  of  analysis  and  report¬ 
writing,  it  must  be  reserved  for  further  development  in  the  Phase 
II  research. 

Now,  let  us  return  to  the  intriguing  SN  path.  This  path  is 
the  only  thing  which  distinguishes  the  cognitive  system  from  a 
mere  cork  on  the  currents  of  the  universe.  There  is  a  technique 
in  the  theory  of  the  Monte  Carlo  method  which  is  called  "impor¬ 
tance  sampling",  in  which  the  experimenter  salts  the  random  data 
with  certain  rare  events  which  are  known  to  have  an  important 
effect  on  the  simulations,  but  which  are  too  rare  to  just  sit  and 
wait  for  in  truly  random  data.  In  a  similar  fashion,  a  cognitive 
system  requires  repetition  in  order  to  ferret  out  the  associa¬ 
tions  and  the  invar iajnts  which  it  needs  to  build  its  internal 
models.  But  novel  events,  by  definition,  do  not  present  them¬ 
selves  at  frequent  intervals.  Therefore,  the  cognitive  system 
must  have  a  way  to  salt  its  observations,  and  it  does  this  by 
manipulating  its  environment  to  repeat  the  novel  event  or  to 
inspect  it  from  a  different  angle  so  that  the  tentative  learning 
(called  a  hypothesis)  can  be  tested  eund  adjusted  before  it  evapo¬ 
rates.  This  is  a  necessary  function  of  the  SN  path.  Essential¬ 
ly,  it  exists  as  a  means  to  "salt"  the  experience  of  the  network 
and  improve  the  efficacy  of  learning.  But  it  can  also  serve  to 
drive  B'  toward  B"  in  the  event  that  learning  fails  to  drive  B" 
toward  B' . 

It  is  tempting  to  suggest  that  the  entropic  structure 
analysis,  via  DMORPH,  can  once  again  be  brought  to  bear  by  using 
it  to  compare  the  structural  similarity  between  the  two  represen¬ 
tations,  B’  and  B",  but  this  would  not  be  appropriate.  DMORPH 
only  applies  to  random  vectors  and  ergodic  stochastic  processes, 
whereas  the  comparison  between  B’  and  B"  which  is  used  for 
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learning  mus'b  be  almosb  instantaneously  computed  on  a  sample-by- 
sample  basis.  It  now  seems  that  the  two  structures  between  which 
DMORPH  might  be  expected  to  find  similarity  are  the  state  tran¬ 
sition  operator  of  the  observed  system  and  the  state  transition 
operator  of  the  neural  system.  In  the  case  of  an  artificial 
neural  network,  this  operator  is  available  in  the  form  of  the 
matrix  of  synaptic  weights  and  the  associated  nonlinear  transfer 
characteristics,  but  they  must  be  treated  as  random  operators 
(cf.,  Skorohod  [81,  for  the  linear  case)  or  else  the  entropy  will 
collapse  to  zero.  If  this  seems  difficult  to  ct.rry  out,  it  is  no 
doubt  far  easier  than  the  other  half  of  the  problem,  which  is  to 
estimate  the  transition  operator  of  the  observed  system  (except 
in  the  most  simple  and  controlled  experimental  conf igurations) . 

It  may  well  be,  as  we  have  alluded  to  above,  that  in  order  to 
measure  the  structural  similarities  in  such  complex  systems,  our 
measurement  instrument  will  have  to  contain  the  cognitive  equiva¬ 
lent  of  the  system  it  is  trying  to  measure. 

In  closing  this  analysis,  we  would  be  remiss  if  we  did  not 
acknowledge  the  role  of  experiment  in  the  development  of  our  own 
cognitive  models  of  cognitive  systems.  We  have  developed  a  tool, 
the  DMORPH  algorithm,  which  provided  corrective  inputs  to  our 
models  both  by  its  success  in  measuring  structural  similarity  in 
random  signals,  and  by  its  success  in  showing  that  such  structure 
is  not  relevant  in  the  ways  that  we  had  presupposed  when  we 
initiated  this  project.  The  paths  to  understanding  are  often 
more  sinuous  than  direct.  We  had  intended  that  DMORPH  should  be 
applied  to  the  reverse-engineering  of  the  brain  in  a  direct, 
gradient  ascent  assault.  Instead,  it  is  in  the  design  and 
verification  of  DMORPH  that  we  found  unexpected  clues  to  the 
organization  of  cognitive  systems.  We  expect  that  this  sinuous 
path  will  now  be  more  productive  than  the  original  plan. 
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I 

8.  CONCLUSIONS  AND  APPLICATIONS 

We  are  confident  that  cognition  is  simply  not  to  be  under¬ 
stood  in  isolation  from  the  essential  interaction  of  the  cogni¬ 
tive  system  with  the  environment  which  it  learns  to  comprehend. 

No  neural  network,  however  complex,  will  exhibit  cognition  if  it 
is  relegated  to  passive  observation  of  its  environment.  The 
conclusion  here  is  perhaps  the  only  solid  confirmation  of  a  pre¬ 
conceived  idea  which  we  had  prior  to  beginning  this  work:  It  is 
that  significemt  progress  with  neural  networks  cannot  be  expected 
without  the  maintenance  of  close  ties  with  biology. 

This  work  shows  the  potential  value  of  stochastic  structure 
analysis  in  the  design  and  improvement  of  neural  network  models 
and  it  is  clear  that  in  six  months  we  have  only  begun  the  process 
of  testing  and  analysis  of  the  various  network  designs.  Now  that 
the  software  tools  are  available,  the  structural  analysis 
deserves  to  be  carried  out  in  a  thorough  and  organized  fashion  on 
many  of  the  existing  network  architectures  to  determine  whether 
and  to  what  degree  their  learning  algorithms  record  and  reproduce 
the  structure  in  the  signals  that  they  observe. 

Although  the  DMORPH  algorithm  meiy  not  be  applicable  as 

originally  planned  to  an  automatic  architecture  mapping  scheme, 

it  clearly  has  an  immediate  utility  for  the  evaluation  and 

improvement  of  neural  network  learning  algorithms  and  transfer 

equations.  Our  intention  is  to  refine  the  algorithm  (its  sort 

routine  could  be  made  faster  and  less  sensitive  to  multiplicity 

of  data)  and  commercialize  it  as  a  utility  to  our  commercial 

TM 

neural  network  simulation  package,  SYSPRO  Furthermore,  it  is 

clear  that  we  can  employ  both  SYSPRO  and  DMORPH  for  neural 
network  design  and  evaluation  for  the  benefit  of  the  Air  Force 
and  other  government  agencies  who  must  compare  designs  and 
determine  their  applicability  to  their  needs. 
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Fig.  A13  Schematic  of  recording  and  conditioning  areas  on  MMEP .  PgAll 

Fig.  A 14  4(X)  neuron  culture  on  recording  area  of  MMEP  1 .  Pg  A12 

Fig.  A  IS  Low  density  culture  on  ITO .  Pg  A12 , 

5.  REPRESENTATIVE  ANALOG  DATA .  Pg  A 13 

Fig.  A16  Multitrace  oscilloscope  representation  of  multichannel  data.  PgA14 

Fig.  A17  Typical  action  potentials .  PgAlS 

Fig.  A18  Olfactory  bulb  activity .  PgA16 

6.  DIGITAL  PROCESSING  SYSTEM .  PgA17 

Fig.  A 19  Present  equipment  setup .  Pg  A 1 8 

Fig.  A20  Present  configuration  of  the  Masscomp  S700  system .  Pg  A19 

Present  Computer  Hardware .  Pg  A20-23 

Data  acquisition  protocols .  Pg  A24 

Fig.  A20  Real  time  display  programs .  Pg  A2S 
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Fig.  A  1.  MMEP  1.  The  lower  picture  shows 
the  entire  5  x  5  cm  plate  with  conductors  and 
amplifier  contacts.  The  center  area  is  enlarged  and 
shown  on  top.  The  36  electrode  arc  arranged  in 
6  rows  and  6  columns  with  200um  and  lOOum  spacing 
respectively.  Each  conductor  is  10  um  wide  and  lum 
thick.  These  plates  were  fabricated  in  1976  free  of 
charge  by  the  Siemens  Corporation  in  Munich,  Germany. 


1'IG.AZ  .Muliiinicrociccirodc  plate  with  64  indium- 
tiii  oxide  electrodes,  (A)  5ctn  x  5cm  x  1.2  nun 
plate  showing  ot'crall  arrangement  of  conductors. 
(B)  Central  region  of  MMI:P  showing  0.5  mm  x  0.5 
mm  recording  matrix  and  radial  conductors.  (C) 
Phase  contrast  microgr.iph  of  the  recording  area 
(conductor  width;  10  pm;  column  spacing:  40  pnt; 
row  spacing;  200;im). 
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Fig.A3Arrangemcnt  of  HO  conductors  in  recording  area  of  inultiniicroelectrode  plate  (MMEP).  The 
plate  was  insulated  with  2  fim  tliick  layer  of  a  polysiloxanc  resin  and  deinsulaled  with  a  single  laser  shot 
at  the  end  of  each  conduettir.  (Jruiind  electrodes  hast*  received  multiple  deinsulation  shots.  Interelectrode 
spacings  for  rows  2-6  are  100  ^rm  between  columns  and  200  pm  between  rows.  In  the  1  X  1.5  mm  center 
recording  area,  the  conductors  are  10  pm  wjde  and  lf\0  nm  thick. 

FROM:  Cross,  G.W.,  W.  Wen  and  J.  Tin  (l')8.5T  Transp.aicnl  indium-tin  oxide  patterns  for 
extracellular,  multisilc  recording  in  neuronal  cultures.  J.  Ncutosci.  Melli.  JL2.:  243-252. 
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2.  MMEP  IMPEDANCES  and  LASER  DE-INSULATION .  Pg  A4 

Fig.  A4  Shunt  capacitance  as  a  function  of  insulation  thickness  ....  Pg  AS 

Fig.  AS  Signal  transfer  as  a  function  of  shunt  impedance .  Pg  AS 

Fig.  A6  Laser  dcinsulation .  Pg  AS 

Fig.  A7  Normal  and  gold  plated  indium-tin  oxide  (ITO)  impedance  as  a 

function  of  crater  area .  Pg  AS 

Fig.  A8  Neurons  on  transparent  ITO .  Pg  A6 

Fig.  A9  Recording  crater  geometries  (ITO) . Pg  A6 
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INSULATION  lO^YER  THICKNESS  (>im] 


riccirodc  sJuml  c;ipacilance  (curv«)  in  pic<dara(I«  on  left 
ordinate  and  shunt  inipcdaitc?  at  1  kllz  (straight  lines)  in  nicgohnis 
on  the  right  ordinate  as  a  function  of  insulation  layer  thickness. 
Calculations  were  carried  out  for  maximum  and  ininimum  con¬ 
ductor  areas  resulting  from  different  conductor  lengths  situated 
under  (he  saline  pool.  Solid  Imes  represent  extreme  values  when 
a  circular  open  culture  chamber  is  utilized  (27  mm  diameter).  Dotted 
lines  result  from  the  tnaximuin  conductor  area  under  saline  when 
a  25  nrm  x  41  niftt  closed  culture  chamber  is  used  (I  ig.  I).  (Relative 
dielectric  constant  =  4.) 


/?r  Percent  of  signal  {E)  seen  by  electrode  lip  rcacliing  atiiplilier 
as  a  function  of  sliunt  impedance  7.^  and  electrode  impedance  7.^ 
(1  klfz).  Serious  signal  attenuation  occurs  at  electrode  impedance 
above  5  Mil  and  shunt  impedances  below  30  Mil.  Only  small  im¬ 
provements  in  electrode  performance  can  be  expected  from  increasing 
the  shunt  impedance  above  50  Mf7.  Curves  were  calculated  with  an 
amplifier  input  impedance  (Z^)  of  20  Mfl  (open  circles)  and  15  MSI 
(solid  circles,  with  equal  to  5  Mfl). 


Figs.  A4  -  A6  From:  Gross,  G.W.  (1979), 
Simultaneous  single  unit  recording  in  vitro 
with  a  photoetched  laser  deinsulated  gold 
multi-microelecirode  surface.  lEEEE  Trans, 
Biomed.  Eng.  BME-26:  273-279. 


ACp  I..nscr-imlucc()  cIcctroJc  diinsulalion  and  concomitant  imped¬ 
ance  ch.inpe  at  I  kllr..  (A)  Intact  gold  conductor  12  um  wide  and 
2  urn  tliick  covered  witli  a  3-4  urn  thick  layer  of  insulation.  (H) 
After  single  laser  slrol  (337  iim,  1  X  10>®  W/cm7)  removed  insula¬ 
tion  fragments  and  gold  p.irlicles  can  be  seen  in  vicinity  of  electrode 
tip.  (C)  Change  in  magnitude  of  t  ktli;  sinusoidal  signal  across 
electrode  at  moment  of  laser  exposure  (arrow).  (D)  Similar  signal 
displayed  after  half-wave  rccliric.ation  on  cliart  recorder.  Electrode 
impedance  decreases  from  42  MS7  to  1.6  MSI  with  one  laser  shot 
and  rises  slowly  to  2.2  MSI  wiitiin  3.5  rntn. 


Fig.  A7  From;  Gross,  G.W,,  W.Wen  and  J.  Lin  (1985). 
Transparent  indium-tin  oxide  patterns  for  extracellular, 
multisile  recording  in  neuronal  cultures.  J.  Neurosci. 
Meth.  lii;  243-252. 


f\n  Kcuording  crater  impedances  as  inverse  fiinetii>ns  of  exposed  area  in  )im’  fru  1  I  O.  gold,  and  IT  O 
that  was  gold  plu'eil  in  tlie  crater.  The  linear  functions  represent  normalircil  impedances  of  11,30  MS2pm^ 
for  ITO  and  255  Nil2/im'  for  gold.  The  higit  impedance  for  the  1 1  O-eleetrolyte  interface  may  be  partially 
a  result  of  further  oxidation  of  the  metal  oxide  during  laser  dcinsulation.  Note  that  gold  plating  of  this 
interface  lowers  impodanees  to  those  established  for  gold  condueliirs. 
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3.  RECORDING  CHAMBER  DESIGN  and  ASSEMBLY. 

Fig.  A 10  Top  and  side  views  of  chamber . 

Fig.  All  Medium  circulation  system . 

Fig.  A12  Chamber  assembly . 


North  Texas  State  University 

as  of  May  1988:  University  of  North  Texas 


PgA7 
PgA8 
PgA8 
Pg  A9 


/17 


Dr.  Guciitcr  W.  Gross  CENTER  for  NETWORK  NEUROSCIENCE 


North  Texas  Stale  University 

as  of  May  1988:  University  of  North  Texas 


I - SSmni - 1 


FIG.  AlOa.  Top  view  of  as.scmblcd  closed  chamber 
and  ciccirodc  plate  holder  showing  the 
mullielectrode  plate  with  its  amplifier  contact 
strips  (Zebra  strips)  to  either  side  and  a  chamber 
cover  containing  the  observation  window. 

Medium  changes  arc  carried  out  via  the  two  ports 
adjacent  to  the  20mm  window. 


FIG.  AlOb.  Side  view  of  closed  chamber  containing 
a  20mm  quartz  or  glass  window  matched  to  the 
objective  to  be  used.  This  arrangement  allows 
la.scr  cell  surgery  with  Zeiss  Ulirafluar  x32  and 
xlOO  objectives  that  have  working  distances  of 
0.45mm  and  0.12mm,  respectively. 


FROM:  Gross.  G.W,  and  M.  Highlowcr  (198(5)  An  approach  to  the  dcicrminaiion  of  network 
propcnics  in  mammalian  neuronal  monolayer  cultures.  In:  Proceedings  of  the  First  IEEE 
Conference  on  Synthetic  Microstruciurcs  in  Biological  Research.  Pcckcrar.  M.C..  Sh.imma. 
S.A.,  and  Wyatt,  R.J.  (cds).  Pp.  3-21.  Washington,  D.C.:  Naval  Rc.search  Uboratoty. 


CLOSED  CIRCULATION  SYSTEM  TO  IMPROVE 
ELECTROPMYSIOLOGICAL  CULTURE  STABILITY 


moist  lOU  CD^ 
tn  oir 


CULTURE  AND  RECORDING 
CHAMBER 

FIG.  A1 1 .  Schematic  drawing  showing  the  recording  chamber  and  closed  circulatory 
system.  The  recording  chamber  contains  about  SOOpl  of  conditioned  medium.  To  maintain 
pH  and  osmolarity  the  medium  in  the  recording  chamber  is  conslanlly  circulated  through 
a  10  ml  reservoir  of  conditioned  medium.  Moi.st  10%  CO2  in  air  is  pumped  into  die  reservoir 

to  maintain  pH.  An  in-line  0.22  pm  filler  in.surcs  sterility  in  the  recording  chamber. 
Phamiacological  agents  can  be  added  to  the  reservoir  and  pumped  to  the  chamber. 
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4.  MINICULTURES  on  MMEPs .  PgAlO 

Fig.  A13  Schematic  of  recording  and  conditioning  areas  on  MMEP .  PgAlO 

Fig.  A 14  400  neuron  culture  on  recording  area  of  MMEP  1 .  Pg  A 1 1 

Fig.  A15  Low  density  culture  on  ITO .  PgA12 
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FIG.AI3  ■  CONFINEMENT  OF  NETWORKS  OVER  RECORDING  AREA  VIA  SELECTIVE 
'  ADHESION. 

The  recording  area  (RA)  is  a  O.^^im  x  1mm  region  in  the  center  of  the 
glass  electrode  date  where  all  conductors  terminate.  Cultures  ^'ist 
be  confined  to  this  area  to  simplify  the  network  analysis.  The  hydro- 
phobic  insulation  material  is  flamed  tlirough  masks  to  generate  specific 
adhesion  patterns.  The  pattern  shovm  consists  of  two  ''conditioi^ng 
areas"  to  either  side  of  a  small  (2mm  diam)  adhesion  island  (Alj 
tered  on  the  recording  area.  The  conditioning  areas  are  necessary  for 
the  proDcr  development  of  neurons.  Neuronal  connections  between  tne 
throe  areas  do  not  develop.  Medium  continuity  exists  at  all  times. 
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Left;  Eig.  A14.  Monolayer  cultures  on  MMEPs.  A 
2  mil)  diainctcr  monolayer  culture  centered  on  die 
recording  matrix  of  a  MMEP  1.  Adhesion  islands  arc 
generated  on  die  nomiully  liydiophobic  insulation 
layer  widi  a  flaming  technique  through  masks.  Laser 
dcinsulation  craters  are  revealed  by  die  habs  at  the 
ends  of  tlic  conductors.  Culture  density:  4(X) 

ticurons/mni2. 

FROM;  Sci.  Aincr.  256:  62. 

Bottom:  Fig.  A  IS.  Center  region  of  a  culture  on  an 
ITO  MMEP  2.  Note  drat  die  transparent  conductors 
do  not  interfere  with  microscopy.  The  heavy  metal 
plating  in  die  crater  is  an  artifact  of  die  Bodian 
histology  method  due  to  prccipiuidoit  of  silver  and 
gold  onto  exposed  ITO.  All  conductors  are  10  pm 
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5.  REPRESENTATIVE  ANALOG  DATA . 

Fig.  A 16  MuUitrace  oscilloscope  representation  of  multichannel  data . 
Fig.  A17  Typical  action  potentials 
Fig.  A18  Typical  burst  activity 
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n *(0  Acliviiy  pailcms  from  one  culture  sampled  34  d  after  seeding.  The  data  show  sequential  oscilloscope  sweeps  (2  sec/div)  stored  at  40 
see  intervals  and  do  not  represent  simultaneous  activity.  In  this  example,  bursting  is  evident  on  at  least  16  of  24  funcliona)  electrodes  or  on  16  of 
20  electrodes  caro  'ng  discernible  activity.  The  nnmhcrs  in  the  n^ht  of  each  panel  rcpicsent  the  row-column  identifiers  of  each  electrode.  E  63 
(sixth  row,  third  column)  and  E  64  were  left  insulated  to  allow  a  monitoring  of  shunt-impedance  stability.  The  symbol  ^  denotes  electrodes  inactive 
because  of  conductor  disconhnuutcs.  //.  Samples  of  complex  smglc-unit  bursting  recorded  from  E  31  over  a  5  min  period  (sweep,  5.0  sec). 

FROM;  Droge.  M.H.,  G.W,  Gross.  M.H.  Hightower,  and  L.E.  Czisny  (1986)  Multiclccirode  analysis 
of  coordinated,  rhythmic  bursting  in  cultured  CNS  monolayer  networks.  J.  Ncurosci.  6: 

1583-1592. 
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I'IG.n  f  I  Characteristic  action  potentials  and  burst 
patterns.  (A)  High  frctnicncy  (SOOlIz)  burst 
showing  decreasing  spike  amplitude.  (U)  Large 
amplitude  (l.GmV)  action  potential  from  the  same 
neuron.  (C)  Small  (140^iV)  spike  rising  from  a 
30)iV  noise  level.  (D&E)  Simultaneous  burst 
patterns  on  two  electrodes.  Note  constant  spike 
amplitude  at  low,  and  decreased  amplitudes  at 
high  firing  frequencies.  Low  amplitude  ionic 
activity  from  a  separate  unit  is  maintained 
between  bursts  in  (L)  Positive  deflection  is  up  in 
all  traces. 


Gross,  G.W.  and  M.L.  Hightower  (1987)  Muliiclectrodc 
investigations  of  network  properties  in  Neural  monolayer 
cultures.  In:  Biomedical  Engineering  ,  Recent  Advances, 
(R.C.  Ebcrhardt,  Ed.),  McGregor  and  Werner,  Washington; 
in  press. 
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Oscilloscope  tracings  from  olfactory  co-culiurcs  (dissociated  olfactory  bulb  cells 
cultured  with  explains  of  olfactory  ncurocpithelium)  grown  on  MMEPs.  Upper  left:  sustained, 
rhythmic  burst  activity;  each  line  contains  ten  seconds  of  data  and  the  average  amplitude  of 
each  burst  is  200p.V  (peak  to  peak).  Upper  right:  single  unit  action  potential;  the  entire 
tracing  represents  20  msec  and  the  amplitude  of  the  action  potential  is  about  ImV  (p/p). 

Lower  left;  single  unit  action  potential  showing  amplitude  decay;  the  entire  tracing 
represents  10  msec  and  the  amplitude  of  tlic  largest  action  potentials  arc  about  800pV  (p/p). 
Lower  right:  expanded  tracing  of  similar  tracings  to  the  right  illustrating  amplitude  decay 
and  waveform  alteration;  the  entire  tracing  represents  50  msec  and  tJic  amplitude  of  the 
largest  action  potentials  arc  about  800p.V  (p/p). 
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6.  DIGITAL  PROCESSING  SYSTEM .  Pg  A17 

Fig.  A19  Present  equipment  setup .  Pg  A18 

Fig.  A20(a)  Present  configuration  of  the  Masscomp  5700  system .  Pg  A19 

Present  Computer  Hardware .  Pg  A20-23 

Data  acquisition  protocols .  Pg  A24 

Fig.  A20(b)  Real  time  display  programs .  Pg  A25 

Examples  of  data  processing .  Pg  A26-29 
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Fig.  A  19  Present  data  acquisition  and  processing  setup.  The  first  stage  amplifiers  reside  on  the 
microscope  stage  to  either  side  to  the  MMEP.  Second  stage  amps  are  connected  to  a  patch  panel,  to  an  LED 
display,  to  an  auditory  monitor  and  to  the  Masscomp  5700  computer.  Integrators,  an  8  channel  strip 
chart,  and  oscilloscopes  arc  serviced  by  the  patch  panel.  The  LED  display  represents  the  physical 
layout  of  the  electrodes  on  the  MMEP  and  displays  activity  in  a  three  color  sequence  (green  -  yellow  - 
red)  depending  on  spike  intensity.  The  auditory  monitor  is  fed  by  the  LED  circuit  and  presents  activity 
on  each  electrode  as  a  different  carrier  frequency.  Teh  later  two  analog  devices  are  very  useful  for 
determining  which  electrodes  arc  active  and  also  for  the  recogniton  of  pailcnis. 


MASSCOMP  5700:  Current  Configuration 
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PRESENT  COMPUTER  HARDWARE 


1.0  SCOPE 

This  document  outlines  the  computational  facilities  available  at  the  Center  for  Network 
Neuroscience.  The  principle  system  is  a  Masscomp  5700  parallel  processing  computer. 

The  Masscomp  5700  is  intended  for  algorithm  development,  statistical  manipulation,  and 
real  time  experiment  monitoring.  The  center  also  has  a  MicroVax  f-r  experiment  control 
and  a  network  of  Macintosh  Computers  for  document  preparation.  i  . 

2.0  MASSCOMP  5700 

2. 1  System  Definition 

The  Masscomp  5700  is  a  computer  mainframe  capable  of  clustering  four  different  types  of 
specialized  processors.  The  mainframe  is  developed  on  the  Motorola  68000  family  of 
computer  processors.  The  processors  include  a  standard  CPU,  a  Data  Acquisition  and 
Control  Processor  (DACP),  a  Pipeline  Processor,  and  a  Graphics  Processor.  The  current 
system  configuration  employs  a  single  DACP  and  two  of  each  of  the  other  processors. 

The  system  provides  access  to  two  industry  standard  busses  for  peripherals  as  well  as  a 
high-speed  main  bus.  Mtiltibus  provides  access  to  standard  computer  peripherals  such  as 
disk  drives  and  tapes.  STD  bus  is  used  for  experimental  instrumentation. 

2.1.1  Functional  Description 

2  1.1.1  The  Standard  CPU 

Masscomp  has  two  different  modules  tJiat  can  be  used  for  the  standard  CPU.  The 
Center  has  two  of  the  68020  modules.  The  68020  module  contains  a  68020  CPU, 
68881  Math  Coprocessor,  an  8K  Cache  area,  and  a  Multibus  Adapter.  The  math 
expansion  module,  the  lightning  floating  point  module,  expands  the 
throughput  of  the  math  co-processor  on  scientific  functions.  Belli  of  the  68020 
modules  are  equipped  with  the  lighting  boards.  Each  processor  is  capable  of  about 

3  Mflops  a  second. 

2. 1.1.2  The  Data  Acquisition  and  Control  Processor 

The  DACP  is  an  8  MHz  bit-slice  processor  that  is  intended  for  realtime 
operations.  The  DACP  is  located  on  multibus  and  provides  an  adapter  to  STD  bus. 

The  DACP  controls  service  inlcrrupls  from  the  STD  bus  modules  and  load  blocks  of 
read  data  into  main  memory.  The  center  currently  has  four  STD  bus  modules,  a 
clock,  j  1  MHz  A/D,  and  two  sample  and  hold  modules. 

2. 1.1. 3  The  Pipeline  Processor 

The  math  pipeline  processor  u.ses  a  7. 1  MHz  adder  and  multiplier  pipe  to 
provide  a  pcrfomiancc  of  14.2  Mflops  per  second.  The  system  has  an 
instruction  queue  for  DMA  operations  as  well  as  for  math  operation.  DMA  can  be 
pcrfomicd  simultaneously  with  math  operations  in  different  sections  of  the  128KB 
of  memory. 
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2. 1.1. 4  Graphics  Processor 

The  graphics  processor  accepts  high  level  graphic  commands  for  an  Aurora 
Display.  The  Aurora  is  an  1 150  x  910  pixel  display  with  4096  displayaul*.,  »,olors  out 
of  a  16  million  color  palette.  The  graphics  processor  also  controls  the  I/O  from  the 
Aurora  keyboard  and  mouse. 

2.1.2  Peripherals 

2. 1.2.1  Memory  ,  . 

The  system  is  configured  with  8  MB  of  main  memory  and  6  MB  of  graphics 
memoiy  located  on  multibus. 

2.1.2.2  Mass  Storage 

The  system  contains  two  Fujitsu  Eagle  disk  drives  as  the  principle  mass  storage 
device.  Each  Eagle  has  387  MB  of  disk  storage.  The  system  is  also  equipped  with  a 
1/2"  tape  drive  for  doing  backups  and  a  5  1/4”  floppy  disk  drive  for  system 
configuration  and  software  updates. 

2. 1.2.3  I/O 

The  system  has  14  RS-232  Serial  Ports  and  an  Ethernet  connection.  The  system 
uses  two  of  the  RS-232  ports  for  interfacing  to  the  NTSU  broadband  network  for 
terminal  access  and  another  R3-232  connection  to  a  1200  baud  modem.  The  system 
is  also  connected  to  two  dot-matrix  printers,  and  a  VTIOO  that  serves  as  system 
monitor. 

2.2  Pcrfonnance 

Each  standard  CPU  has  a  benchmark  of  3300  Kwhetstoncs.  Each  pipeline  processor  is 
capable  of  14.2  Mflops  sustained  throughput.  The  combined  throughput  of  two  standard 
CPU  modules  and  two  pipeline  processors  is  estimated  at  35  Mflops. 

2.3  Physical  Attributes 

The  system  is  housed  in  two  cabinets.  The  primary  cabinet  houses  the  30  slot  frame  for 
Masscomp  bus  and  multibus  boards,  one  Eagle  hard  disk,  and  two  9  slot  STD  bus  frames. 
The  second  cabinet  contains  a  second  Eagle  and  a  tape  drive. 

2.4  Maintenance  and  Support 

The  system  is  maintained  on  a  service  contract  that  provides  for  replacement  of  defective 
hardware,  as  well  as  software  support  and  routine  system  maintenance. 
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3.0  MICROVAX 

3. 1  System  Definition 

The  MicroVax  was  purchased  to  do  processor  control  of  lab  equipment.  The  control 
process  will  be  done  through  a  CAMAC  interface. 

3.1.1  Functional  Description 

The  system  is  a  single  board  version  of  the  Vax  mainframe  produced  by  Digital  Equipment 
Corporation.  The  system  is  equivalent  to  a  1 1/785  with  a  math  co-processor  >^ith  Slightly 
slower  bus  hardware. 

The  system  is  also  equipped  with  a  graphics  processor  with  a  display  of  1024  x  1024  pixels. 

Each  pixel  can  be  assigned  one  of  16  colors  selected  from  a  4096  color  palette. 

3.1.2  Peripherals 

3. 1.2.1  Memory 

The  MicroVax  is  configured  with  3  MB  of  memoiy. 

3. 1.2.2  Mass  Storage 

The  MicroVax  has  a  70  MB  Winchester  disk  drive  with  a  T50  streaming  tape 
backup. 

3.1.2.3  170 

The  MicroVax  has  two  RS-232  ports  and  an  Ethernet  connection. 

3.2  Performance 

The  MicroVax  is  benchmarked  at  1000  kilowhetstones. 

3.3  Physical  Attributes 

The  MicroVax  fits  into  a  standard  rack  mount  cabinet  using  5"  of  space. 

3.4  Maintenance  and  Support 

The  system  is  being  serviced  as  needed,  with  billings  for  time  and  materials.  Software 
support  is  being  provided  through  campus  computer  center. 
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4.0  MACKINTOSH  NETWORK 

4. 1  System  Elefinition 

The  System  is  a  network  of  six  workstations  connected  to  a  laser  printer.  The  system 
provides  the  Center  with  document  formating  from  personal  computers. 

4.1.1  Functional  Description 

4.1. 1.1  Network  Description 

The  Network  uses  Applenet  to  connect  the  work.stalions  and  laser  printer 
together.  Applenet  is  a  broadband  network  similar  to  ethemet  with  a  1/4  Mhz 
bandwidth. 

4. 1.1. 2  Workstation  Description 

The  center  has  five  Mac+  workstations  and  a  Mac  2  workstation.  The  Mac+  is  a  1MB 
system  with  a  black  and  white  screen.  The  Mac  2  is  a  color  system  with  16  color 
pixels  and  a  40  MB  Winchester  hard  disk. 

4. 1.1. 3  Printers 

The  main  printer  is  an  Apple  LaserWriter  printer.  The  LaserWriter  is  capable  of 
printing  at  the  rate  of  8  pages  a  minute.  Two  workstations  have  2  dot  matrix 
printers  attached  for  rough  drafts. 

4.2  Maintenance  and  Support 

Hardware  failures  are  fixed  as  needed.  Software  support  is  supplied  by  the  vendor. 
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Fig.  \30  Enicigiiig  rcallinic  display  program  for  3d  electrodes  showing  spike  dala  (top)  and  integrated  data 
(bottom).  All  electrodes  arc  displayed  in  the  lower  right  of  the  screen  in  a  physical  arrangement  that  mimics  the 
layout  of  the  electrodes  under  the  microscope  (MEI’  I),  lire  center  of  the  stpiarc  electrode  selection  buttons  will 
show  different  colors  as  a  function  of  burst  intcirsity  (not  com(ilcted).  1  he  squatc  collar  has  a  color  code  that 
corresponds  to  die  trace  position  on  die  screen.  Sweep  speeds  and  amplitudes  arc  selected  at  die  bottom  of  the 
screen.  A  touil  of  6  channels  will  eventually  appear  on  the  screen.  .Sl.alistical  parameters  such  as  average  burst 
duration,  l,..ist  periud,  buist  aiea,  and  Imrst  l\(ic  (as  well  ns  v.iiioiis  hisioj-.i.iiiis)  will  also  be  scleclablc  from 
the  piuicl  alxive  the  cicctrude  nialiis.  Most  paianicteis  can  he  |>loilcd  as  a  luiielion  of  lime  and  will  provide  a 
continuous  6  hour  rceoid  of  llicsc  pai  anielcrs. 

Because  of  our  constant  csposiiie  to  niuhicliaiincl  d.'il;i.  we  li;>vc  a  high  piohahility  of  developing  realistic, 
effective,  opcralor  fticndly  piograins  that  mighi  Ic  n'cd  lor  iii  iiiy  dillcrciil  iiiiillicliaiiiicl  piohlcms.  It  is 
conceivable  that  ECU  dala  could  be  selected  in  die  .s.iiiie  v\ay  Iroiii  a  hiaiii/clccliodc  schematic  that  replaces  die 
6x6  electrode  matrix. 


Dr.  Guenter  W.  Gross 


CENTER  for  NETWORK  NEUROSCIENCE  North  Texas  Slate  University 

as  of  May  1988;  University  of  North  Texas 
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A  second  data  processing  example,  exhibiting  a  liiglicr  level  of  A/C  noise.  After 
filtering  for  line  noise,  die  R/C  integration  of  the  binary  channel  produces  an 
electrically  "clean”  signal. 
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Source 

File; 

DTEST.FOR  Options;  /ELY  02/03/88  16 

1 

PROGRAM  DTEST 

2  C 

3  C 

PROGRAM  TO  TEST  THE  "GIBBS"  SUBROUTINE  (RLD  12/02/87 i 

4  C 

5 

DIMENSION  X(8,2048hY(8,2048),  TILEX(256),  TILEY(2565 

6 

INTEGER*2  HIST0( 0: 255, 0: 255) .  INDEX(2048) 

7  C 

8 

CHARACTER*8  R4NAMES.  I4NAMES,  LINAMES 

9 

REAL*4  R4VALUE 

10 

INTEGER+4  I4VALUE 

11 

L0GICAL>('1  LIVALUE,  QUIT,  PLOT 

12  C 

13 

COMMON  /NLIST/  R4NAMESf 21 ) ,  I4NAMES(21 ) ,  L1NAMES( 7 ) , 

14 

1 

R4VALUE(21) , I4VALUE( 2 1 ) .  L1VALUE(  7 ) 

15  C 

16 

EQUIVALENCE  ( SEED, R4 VALUED  1 ) ) ,  (A, R4VALUE(2) ) ,  (B, R4VALUEf 3 ) ) 

17 

1 

^GOTIME,R4VALUEM)  ) 

IS 

EQUIVALENCE  (NX. I4VALUE(1)), (NY, I4VALUE(2)). (LEXP, I4VALUE(3 

19 

2 

( INSTR. I4VALUE(4) ), ( lOUNIT, I4VALUE(5) ) , ( IC, I4VALUE(6) )  . 

20 

2 

( IPR, I4VALUE(7) ) , ( IDU, I4VALUE(8) ) . ( IFUN. I4VALUE(9) ) , 

21 

2 

(LDBUG, I 4 VALUE ( 10) ) 

22 

EQUIVALENCE  (QUIT. L1VALUE( 1 ) ) 

23  C 

24  C 

25 

NX  =  8 

26 

NY  =  8 

27 

IC  =  8 

28 

LEN  =  255 

29 

SEED  =  2.468E+12 

30 

A  =  0. 

31 

B  =  1. 

32 

LEXP  =  IC  *  2**MAX(NX,NY) 

33 

INSTR  =  1 

34 

QUIT=. FALSE. 

35 

lOUNIT  =  6 

36 

IDU  =  0 

37 

IPR  =  0 

38 

IFUN  =  0 

39 

LDBUG  =  0 

40  C 

41 

R4NAMES(1)  =  'SEED' 

42 

R4NAMES(2)  =  'A' 

43 

R4NAMES(3)  =  'B' 

44 

R4NAMES(4)  =  'GOTIME' 

45  C 

46 

I4NAMES(1)  -  'NX' 

47 

I4NAMES(2)  =  'NY' 

48 

I4NAMES(3)  =  'LEXP' 

49 

I4NAMES(4)  =  'INSTR' 

50 

I4NAMES(5)  =  'INPTUNIT' 

51 

I4NAMES(6)  =  'IC' 

52 

I4NAMES(7)  =  'PRNTUNIT' 

53 

I4NAMES(8)  =  'DATAUNIT' 

54 

I4NAMES(9)  =  'IFUN' 

55 

I4NAMES(10)=  'LDBUG' 

56  C 

57 

LlNAMES(l)  =  'QUIT' 

58  C 

59 

NLl  =  4 

60 

NL2  =  10 
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61 

62 

63  C 

64  10 

65  C 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75  C 

76 

77 

78 

79  11 

80  12 
81 

82  C 

83  C 

84 

85  C 

86 
87 
83 

89 

90 

91 

92  C 

93 

94 

95  20 

96 

97 

98  902 

99 

100  C 

101 
102 

103 

104 

105 

106 

107  C 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117  101 

118 

119 

120  100 


NL3  =  1 
TIME=0. 

CONTINUE 

CALL  NAMELISK  lOUNIT,  NLl .  NL2.  NL3  ) 

IF (QUIT)  STOP 
ASEED  =  SEED 
NHX  =  2**NX-1 
NHY  =  2**NY-1 
NTRIALS  =  0 

IFdDU.GT.O  .AND.  GOTIME.  LT.  TIME)  REWIND(IDU) 

TIME  =  GOTIME 

IF( INSTR. EQ. 2)  THEN 

WRITE(9,  11)  (R4NAMES( I) ,  R4VALUE( I) , 1=1,  NLl ; 

WRITER  9, 12)  (I4NAMES(I). I4VALUE(I), 1  =  1, NL2) 

F0RMAT( (4(A8.  ’=  ’,F8.4,1X))) 

FORMATC (4(A8, ’=  ’.18. IX))) 

END  IF 

GENERATE  OR  READ  THE  DATA 

PRINT  ♦,  ’  DTEST :  Generating  raindom  data  for  X  and  Y.  ’ 

DO  20  ITRIAl,  =  l.LE.XP 
IF(IDU.LE.O)  THEN 

CALL  QRX(X( 1,  ITRIAL) .NX. Y( 1, ITRIAL) , NV, A, B. ASEED. IFUN) 

ELSE 

STOP  ’DTEST:  Real  data  initialization  not  available.’ 

CALL  RDATA(X.N,  TIME, IDU) 

END  IF 

CONTINUE 

IFfLDBUG.  GE.  3)  WRITE{6.g02)  (  { X(  I ,  J )  ,  1  =  1 ,  NX ) ,  J  =  1 ,  LEXJ’ ) 

FORMAT! 8 ( IX, F8. 5) ) 

GENERATE  THE  X  AND  Y  EVENT-SPACE  TILINGS 

PRINT  ♦, ’  DTEST:  Generating  the  X  event-space  tilings. ’ 

CALL  UNIVENT!X.  NX,  TILEX, INDEX, IC.LDBUG) 

PRINT  *, ’  DTEST:  Generating  the  Y  event-space  tilings.’ 

CALL  UNIVENT! Y,  NY,  TILEY, INDEX, IC, LDBUG) 

COMPUTE  THE  NORMALIZED  STRUCTURE  INDEX. 

PRINT  '  DTEST:  Computing  the  normalized  structure  index.’ 

DO  100  ITRIAL  =  l.LEXP 

CALL  GIBBS!X! 1,  ITRIAL) , NX, TILEX,  Y! 1 , ITRIAL) . NY, TILEY,  NTRIALS, 
+  HISTO,  NHX,  NHY,  INSTR.  HX, HY. H, G, LDBUG ) 

IF! INSTR. EQ. 2)  THEN 

WRITE!9, 101)  H.HX.HY.G 
FORMAT!4!F12. 5,  IX) ) 

END  IF 


CONTINUE 
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121  C 

122 

123 

124 

125  C 

126 

127  C 

128 

129  901 

130 

131 

132 

133  C 

134 

135 


NTR  =  NTRIALS 

CALL  GIBBS(X, NX.TILEX,  Y.NY.TILEY,  NTR. 
h  HISTO.  NHX.NHY,  0.  HX, HY. H. G . LDBUG ) 


IF(IPR.NE.6)  WRITERS, 901)  NTRIALS. SEED, NX, NY, H. HX, HY. G 


IF(IPR.GT.O)  WRITE( IPR, 901)  NTRIALS, SEED, NX. NY, H . HX, HY, G 
FORMAT(’  DMORPH  EXPERIMENT:  #  TRIALS  =  ’.15,’,  SEED  =  ’, 

1  F10.7.’.  X-DIM  =  M2.'.  Y-DIM  =  ’,12. 

2  /8X, ’WHOLE  ENTROPY  =  ’.F10.7,’,  X  ENTROPY  =  ’,F10.7, 

3  /8X,’Y  ENTROPY  =  ’.F10.7.’.  DMORPH  =  ’,F10.7./) 


GO  TO  10 
END 


NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


136 

137 

SUBROUTINE  QRX(X.NX.Y.I 

138 

DIMENSION  X(NX),Y(NY) 

139 

C 

140 

DO  100  I  =  l.NX 

141 

X(I)  =  A  +  C 

142 

100 

CONTINUE 

143 

C 

144 

DO  200  I  =  l.NY 

145 

Y(I)  =  A  +  ( 

146 

200 

CONTINUE 

147 

148 

GOTOdOl,  102,  103,  104) 

149 

GOTO  1000 

150 

C 

151 

101 

CONTINUE 

152 

C 

X(4)  =  Y(2)  +  X(4) 

153 

GOTO  1000 

154 

C 

155 

102 

CONTINUE 

156 

X(4)  =  (Y(2)+X(4) )/2 

157 

GOTO  1000 

158 

159 

103 

CONTINUE 

160 

DO  1031  1=1. NX 

161 

1031 

X(I)  =  Yd) 

162 

UOTO  1000 

163 

164 

104 

CONTINUE 

165 

X(2)  =  X(l) 

166 

167 

1000 

RETURN 

168 

END 
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NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT;  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


169 

170 

171 

C 

172 

C 

173 

C 

174 

C 

175 

C 

176 

C 

177 

C 

178 

C 

179 

C 

180 

C 

181 

182 

183 

184 

185 

10 

186 

187 

188 

901 

189 

190 

20 

191 

192 

193 

194 

100 

195 

196 

197 

193 

199 

200 

201 

1100 

202 

203 

1200 

204 

205 

206 

207 

SUBROUTINE  RDATA(X.N.T, IDU, ITRIAL) 

This  subroutine  reads  sample  data  from  the  file  FORTn,  where 
n  =  IDU.  It  skips  all  records  with  time-tags  less  than  T, 
reads  all  records  with  time-tags  equal  to  the  time-value  first 
encountered  which  is  greater  than  or  equal  to  T,  and  puts  the 
next,  value  of  the  tirne-tag  into  the  T  variable  before  returning. 

If  you  want  to  interpolate  or  extend  the  data  between  times 
existing  on  the  fi le.  you  must  do  that  external  to  this  subroutine. 

DIMENSION  X(N) ,M(4) ,R{4) 

IFaTRIAL.GT.  1)  GOTO  20 

CONTINUE 

READ(IDU.901.END=1200)  TT.NX.GX, (M(I).R{I). 1=1,4) 

F0RMAT(F9. 3, IX. II, 1X,F8. 3, IX.  4(13,  1X,F10. 3. IX)) 

CONTINUE 

IF(TT.LT.T)  GOTO  10 
DO  100, 1=1, 4 

IF(MfI).GT.0  .AND.  M(I).LE.N)  X(M(I))=R(I) 

TTl  =  TT 

READ ( IDU, 901, END=1 100)  TT.NX.GX, { M( I ) , R(I ) , 1=1 . 4 ) 

IF(TT.EQ.TTl)  GOTO  20 
T  =  TT 

CONTINUE 

RETURN 

CONTINUE 

PRINT  *,  '  NO  DATA  EXISTS  ON  INPUT  DATA  UNIT  MDU 
PRINT  ’  BEYOND  THE  REQUESTED  TIME  T  =  ',T 
STOP 
END 


NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


NUMBER  OF  WARNINGS  IN  COMPILATION  ;  0 
NUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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1  SUBROUTINE  UNIVENT(X.  NX.  THRESH.  INDEX,  IC,  LDBUG) 

2 

3  C*** ****♦♦*****■**+;*  +  *♦* ♦+*****♦**■**♦************ **********H<*:**^*  + ***»)(<;♦*» 

4  C  This  routine  determines  the  SUM  (i=l,K)  2*t(I-l)  threshholds 

5  C  which  will  divide  the  k-dimensional  space  into  regions  that 

6  C  will  have  an  equal  amounts  of  counts  when  the  sample  is 

7  C  drawn  from  the  same  underlying  distribution  which  generated 

8  C  the  c^2*+K  vectors  used  by  this  program  to  set  the  boundaries. 

9  C  If  IC=4,  K=8  and  the  data  is  found  on  lUNIT  =  1  then  the  num. 

10  C  of  threshholds  which  need  to  be  found  are  1+2+4+8+16+32+64+128 

11  C  or  255  based  on  the  512  data  vectors. 

13 

14  + +  +  +  **+■*+♦***♦**** >1' +=+**  +  ***♦**♦♦****♦**+*******;+*+* * 

15  C 

16  C  IC  is  the  integer  multiple  of  the  min.  of  samples ( 2*^NX) . 

17  C  NX  is  the  sample  vector  dimension. 

18  C  NSAMPLE  =  IC  *  2  **  NX,  is  the  number  of  samples. 

19  C  NTHRESH  =  2  NX  -  1,  is  the  total  number  of  thresholds. 

20  C  DATA(I,J),  i=l,NX  ),  J=l, NSAMPLE  )  is  the  sampled  data. 

21  C  THRESH  is  the  array  in  which  the  thresholds  are  stored. 

22  C  LTHR  is  the  length  of  THRESH  and  must  =  -1+2^*NX 

23  C  INDEX  is  a  workspace  integer  array  of  length  NSAMPLE. 

24  C  , 

25  C  The  calling  sequence  for  UNIVENT  is  as  follows; 

26  C  CALL  UN I VENT (  X,  NX,  THRESH,  LTHR,  IC  ) 

27  C 

^ 

29 

30  DIMENSION  THRESH(*) 

31  INTEGER*2  INDEXC*) 

32 

33  C  NOTE'  The  first  X-dimension  (below)  MUST  be  exactly  the  same  as 

34  C  as  in  the  calling  program,  even  if  NX  may  be  different! 

35 

36  DIMENSION  X(8,*) 

37 

38  C  If  your  data  is  integer,  remove  the  comment  from  column  1 

39  C  of  the  next  line  of  code. 

40  C  INTEGER  DATA,DATAT 

41 

42 

43  NSAMPLE  =  IC  *  2  NX 

44  NTHRESH  =  2  ♦*  NX  -  1 

45  C 

46  C  Initialize  the  index  array. 

47  C 

48  IF(LDBUG.GE. 1)  PRINT  +,  ’NSAMPLE  ’.NSAMPLE,’  NTHRESH  ’.NTHRESH 

49 

50  DO  10  I  =  1  ,  NSAMPLE 

51  INDEX(I)  =  I 

52  10  CONTINUE 

53 

54  c** ♦♦**;t*;t^*’i'*  +  ***^**** *♦**♦**♦  ♦♦♦♦t*****’t‘*****^t*  +  t*’t*^*^;<‘*^>t‘*'**>f^t!t* 

55  C 

56  C  For  each  of  the  NX  dimensions,  I,  of  the  data  vector,  we 

57  C  determine  the  2**(I-1)  thresholds  which  divide  the  space  into 

58  C  approximately  equal  (based  on  the  sample)  probability  bins  given 

59  C  that  we  have  already  divided  the  space  for  all  dimensions  less 

60  C  than  I  and  we  consider  for  each  of  the  2**(I-1)  thresholds 
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61  C  only  the  sample  data  in  one  of  the  previous  2t'f'(I-l-l)  bins. 

62  C  For  example,  if  1=1,  all  of  the  sample  is  divided  into  one  of  two 

63  C  bins  based  on  the  value  of  the  1st  component  of  the  NSAMPLE 

64  C  N.X-vectors.  Based  on  this  2**fl-l)=l  threshold  and  the 

65  C  value  of  the  2nd  component  of  the  data  vectors  in  each  of  the  two 

66  C  subsets  of  the  NSAMPLE  data  points,  these  subsets  are  then  divided 

67  C  with  2**(2-l)=2  thresholds.  The  threshold  used  is  the  median  of 

68  C  the  sampled  component. 

69  C 


71 

72 

NS  =  1 

73 

IC2PNXI  =  IC 

*  2**!NX+1) 

74 

DO  130  I  =  1  , 

NX 

75 

76 

JSTOP  =  0 

77 

IC2PNXI  = 

IC2PNXI/2 

78 

JSTART=  1- 

IC2PNXI 

79 

30 

81 

DO  120  L  = 

1  .  NS 

82 

83 

JSTART 

=  JSTART  +  IC2PNXI 

84 

JSTOP 

=  JSTOP  +  IC2PNXI 

85 

86 

87 

IF!LDBUa.GE. 3) 

THEN 

88 

WRITE 

!6,t>  'JSTART  ’.JSTART,’  JSTOP 

89 

WRITE 

! 6, 1002 )  ! INDEX! M) ,  M=1 , NSAMPLE 

90 

END  IF 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 
110 

119 

120 


C  This  code  sorts  the  data  array  IDATA  w.r.t.  it's  Ith 

C  component  (column)  for  the  data  values  corresponding  to 

C  INDEX(JSTART) . INDEX( JSTART+1) .  ....  INDEX( JSTOP) . 

DO  110  J  =  JSTART+1,  JSTOP 

INDEXT  =  INDEX  (J) 

DAT AT  =  X(I, INDEX! J) ) 

IFCLDBUG.GE. 4) 

1  WRITE  (6.*)  'J  ’,J,'  INDEXT  ’, INDEXT, ’  DAT AT  ' , DATAT 

DO  100  K  =  J-l,  JSTART,  -1 
IFCLDBUG.GE. 5)  THEN 

WRITE  (6,*)  'K,  DATAT,  INDEX(K)  K, DATAT, INDEX(K) 
WRITE  (6,*)  'Xd.  INDEX(K) )  ’  .  X(  I ,  INDEX(K) ) 

END  IF 

IF  (  DATAT  .LT.  XC I. INOEX(K) )  )  THEM 

INDEX(K+1)  =  INDEXdO 
INDEX! K)  =  INDEXT 

IF(LDBUG.GE. 5)  THEN 

WRITE  (6,*)  'TEST  DATA  <  DATA  ABOVE, BUBBLE  UP’ 
WRITE  !6,*)  'J  ’,J,’  INDEX! J)  ’, INDEX! J) 
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121  WRITE  (6.*)  'K  ',K, ’  INDEX(K)  ’,INDEX(K) 

122  END  IF 

1  P'1 

124  ELSE 

125 

126  IFfLDBUG.GE. 5)  WRITE  (6.*)  'TEST  DATA  >  DATA  ABOVE,  NEXT  J' 

127 

128  GOTO  110 

129 

130  END  IF 

131 

132  100  CONTINUE 

133 

134  110  CONTINUE 

135 

136  C***^***^*  +  =<'+*  +  **-+^+^*  +  + 1 +  ^*+*+*'t<'******'^'*  +  *  +  **+-+***t^*'^****>f‘t-+******>t'*** 

137  C  Now  use  the  sorted  data  to  determine  the  median/threshold. 


138  C* *♦**+♦  + *+*****=<' ♦'♦■*^*'+***^*****=*'+**  +  *  +  **+-**+*+^** ♦****  +  +  ***  +  **■* ****♦*>!■* 


139 

140 

141 

142  C 

143  C 

144  C 

145  C 

146  C 

147  C 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 


MIDINDEX  =  INT  (  (  JSTOP  +  JSTART  )  /  2  ) 

If  your  data  is  integer,  remove  the  comment  from  column  1 
of  the  next  three  lines  of  code,  and  comment  out  the 
non-floated  definition  of  THRESH. 

THRESH! 2*+(I-l)+L-ll= 

+  (  FLOAT  (  X(I.  INDEX(MIDINDEX)  )  )  + 

+  FLOAT  (  X(I.  INDEX! MIDINDEX+1)  }))/  2. 

THRESH!2tt!I-l)+L-l)= 

+  !  X!  I.  INDEX!MIDINDEX)  )  + 

+  X!  I,  INDEX!MIDINDEX+1)  ))/  2. 

IF!LDBUG.GE. 3)  THEN 

WRITE  !6.*)  'THRESH! 2**! I-l )+L-l,  , 

+  THRESH!2**! I-l)+L-l) 

WRITE  !6.*)  ’  MIDINDEX  ’.MIDINDEX 

WRITE  !6,1002)  !INDEX!M),  M=1,NSAMPLE  ) , I, 2**! I-l )+L-l 
WRITE  ! 6, 1001)  ! THRESH! M),  M=1.NTHRESH  ) 

END  IF 

120  CONTINUE 
NS  =  2tNS 

130  CONTINUE 

IF!LDBUG.GE. 2)  THEN 

WRITE  !6.1001)  !THRESH!M),  M=1,NTHRESH  ) 

WRITE  !6.1004)  ! ! I . J, X! J,  INDEX! I ) ) .  •T=1,NX  ) , 1  =  1 . NSAMPLE  ) 

END  IF 


1001  FORMAT! 10! IX. F7. 3)) 

1002  F0RMAT!20!1X,  13)) 

1003  FORMAT!  F9.0,  IIX,  4  !  IX,  13,  IX,  FIO. 4  )  ) 

1004  F0RMAT!5!1X. 13,  IX, T3.  IX,  F7. 4)) 

RETURN 

END 
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NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


NUMBER  OF  WARNINGS  IN  COMPILATION  ;  0 
NUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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1 

2 

3 

4  C 

5  C 

6  C 

7  C 

8  C 

9  C 

10  C 

11  C 

12 

13 

14 

15 

16 

17 

18 

19 

20  C 

21 
22 

23 

24 

25 

26 

27 

28  10 

29 

30  C 

31 

32 

33 

34 

35 

36 

37 

38  20 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50  C 

51 

52 

53 

54 

55 

56  310 

57 

58  C 

59 

60 


SUBROUTINE  GIBBSf X. NX, TILEX, Y. NY, TILEY, NTRIALS, 

HISTO,  NHX,  NHY,  INSTR,  HX,  HY,  H,  G,  LDBUG) 


This  subroutine  computes  three  entropies  associated  with  two 
random  vectors  X  and  Y,  of  dimensions  NX  and  NY.  H  is  the 
entropy  of  the  concatenated  vectors  after  NTRIALS  of  the  ex¬ 
periment.  HX  and  HY  are  the  separate  entropies  of  X  and  Y 
after  NTRIALS.  G  =  HX+HY-H  is  the  Gibbs  relative  entropy 
of  the  combined  system.  All  entropies  are  computed  with  res¬ 
pect  to  the  tiling  of  the  event  space  specified  by  the  TILE 
arrays . 

DIMENSION  X(NX) , Y(NY) , TILEXC*) , TILEY(*) 

INTEGER*2  HISTO( 0 : NHX, 0: NHY) 

PARAMETER  (ALN2=0. 6931471 ) 

IF{ INSTR  .EQ.  0)  GO  TO  310 
NTRIALS  =  NTRIALS+1 

IDENTIFY  THE  X-EVENT  NUMBER 

KX  =0 
LEVEL  =  1 

DO  10  J=1,NX 

IF(  XfJ)  .GT.  TILEX (LEVEL+KX)  )  KX  =  LEVEL  +  KX 
LEVEL  =  2 ♦LEVEL 
CONTINUE 

IDENTIFY  THE  Y-EVENT  NUMBER 

KY  =0 
LEVEL  =  1 

DO  20  J=1,NY 

IFf  Y(J)  .GT.  TILEY(LEVEL+KY)  )  KY  =  LEVEL  +  KY 
LEVEL  =  2tLEVEL 
CONTINUE 

IF(LDBUG.GE. 1)  THEN 
PRINT  ♦, 'X  =  ’ 

PRINT  ♦, (X(I), 1=1, NX) 

PRINT  ♦, ’  X-EVENT  I.D.  =  ’  ,  KX 
PRINT  ♦ 

PRINT  ♦, 'Y  =  ' 

PRINT  ♦, (Y(I), I=1,NY) 

PRINT  ♦, '  Y-EVENT  I.D.  =  ’ . KY 
END  IF 

BUMP  THE  COUNT  FOR  THE  IDENTIFIED  COMPOSITE  EVENT 
HISTOrKX.KY)  =  HISTO(KX.KY)  +  1 
IF(  INSTR  .EO.  1  )  RETURN 
CONTINUE 

COMPUTE  THE  ENTROPIES  ASSOCIATED  WITH  THE  ACCUMULATED  HISTOGRAM. 
H  =  0. 
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61 

HX=  0. 

62 

HY=  0. 

63 

FLN  -  FLOATfN TRIALS) 

64 

65 

DO  330  1=0. NHX 

66 

MARX  =  0 

67 

DO  320  J=0,NHY 

68 

P  =  FLOAT(HISTO(I,J))/FLN 

69 

IF(P.GT.O)  H  =  H  -  P*LOG(P) 

70 

MARX  =  MARX+HISTO( I, J) 

71 

320 

CONTINUE 

72 

PX  =  FLOAT (MARX) /FLN 

73 

IF(PX.GT.O)  HX  =  HX  -  PX*LOG(PX) 

74 

330 

CONTINUE 

75 

76 

DO  350  J=0,NHY 

77 

MARY  =  0 

78 

DO  340  1=0, NHX 

79 

MARY  =  MARY+HISTOd,  J) 

80 

340 

CONTINUE 

81 

PY  =  FLOAT (MARY) /FLN 

82 

IF(PY.GT.O)  HY  =  HY  -  PY^LOG(PY) 

83 

350 

CONTINUE 

84 

, 

85 

C 

COMPUTE  THE  NORMALIZED  STRUCTURE  FUNCTION 

86 

0 

G  =  2.*(HK+HY  -  H>/( {NX+NY)*ALN2)  [ 

SUPERCEDED  3 

87 

HX  =  HX/ALN2 

88 

HY  =  HY/ALN2 

89 

H  =  H  /ALN2 

90 

SUPHXHY  =  FLOAT (NX+NY) 

91 

AMINH  =  AMAXKHX.HY) 

92 

G  =  (HX+HY-H)/(SUPHXHY-AMINH) 

93 

94 

IF(  INSTR  .EQ,  2  )  RETURN 

95 

96 

C 

RESET  THE  HISTOGRAM  TO  ZERO  FOR  THE 

NEXT  EXPERIMENT 

97 

98 

DO  420  1=0, NHX 

99 

DO  410  J=0,NHY 

100 

HISTOd.J)  =  0 

101 

410 

CONTINUE 

102 

420 

CONTINUE 

103 

NTRIALS  =  0 

104 

105 

RETURN 

106 

END 

NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


NUMBER  OF  WARNINGS  IN  COMPILATION  :  0 
NUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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SUBROUTINE  NAMELIST ( lOUNIT, Nl, N2, N3) 


The  following  declarations  are  for  local  variables. 
CHARACTER  LINE*72, NAME*8. NVAL*72 


9  C 

10  C 

11  C 

12 

13 

14  C 


This  program  is  intended  to  approximate  tho-  NAMELIST  capability 
which  some  FORTRAN  compilers  have,  but  which  RM-FORTRAN  does  not 
have.  The  calling  program  needs  to  have  a  common  block  labeled 
/NLIST/  and  containing  six  arrays: 

CHARACTER*8  R4NAMES 
REAL*4  R4VALUE 


CHARACTER*8 

INTEGER*4 


I4NAMES 

I4VALUE 


18 

19 

20  C 

21 
22 

23  C 

24 

25 

26 

27 

28 

29 

30  C 

31  1 

32  C 

33 

34 

35 

36 

37 

38 

39 

40 

41  2 

42  3 

43  4 

44 

45 

46 

47 

48 

49 

50  C 

51  10 

52  C 

53 

54  12 

55 

56 

57 

58 

59  C 

60  C 


CHARACTER*8 

LOGICAL*! 


LINAMES 
LI VALUE 


COMMON  /NLIST/  R4NAMES( 21 ) . I4NAMES( 21 ) . L1NAMES( 7 ) , 

R4 VALUE ( 2 1 > , I 4VALUE ( 2 1 ) , L 1 VALUE ( 7 ) 

IF{N1.GT. 21. OR. N2.GT. 21 . OR. N3. GT. 7 )  THEN 

PRINT  *, 'Nl  or  N2  or  N3  is  too  large  for  NAMELIST. ' 

PRINT  ♦.’Increase  dimensions  in  /NLIST/  common,  and' 
print  ♦, ’increase  limits  in  first  statement  of  NAMELIST.’ 

STOP 
END  IF 

CONTINUE 

IF(I0UNIT.EQ.6)  THEN 

PRINT  *, ’ENTER  VARIABLE  NAMES  FOLLOWED  BY  VALUES  ACCORDING  TO’ 
PRINT  *,’THE  SYNTAX,  name  =  value  <CR>.  SPACES  ARE  OPTIONAL.’ 
PRINT  *, 'LEGAL  NAMES  AND  CURRENT  VALUES  ARE:’ 

PRINT  * 

WRITE(6,2)  (R4NAMES(J),R4VALUE(J), J=1.N1) 

WRITE(6,3)  (I4NAMES(J>, I4VALUE{J), J=1,N2) 

WRITE(6,4)  (L1NAMES(J),L1VALUE(J), J=1,N3) 

FORMATS (4(2X,A8,  ’ C ’ , F8. 3, ’ 3 ’ ) ) ) 

FORMATS (4(2X,A8,  ' [  ’,16,’  ]’>)) 

FORMAT((6{2X,  A8. ’[  ’,L1,’  ]’)>) 

PRINT  * 

PRINT  *,'IF  YOU  GOOF,  JUST  RE-ENTER  THE  LINE  CORRECTLY.’ 

PRINT  ♦,’ANY  LINE  NOT  HAVING  THE  =  SIGN  IN  IT  TERMINATES  ENTRY,’ 
PRINT  *, 'EXCEPT  •■?"  DISPLAYS  VALUES  AND  STOPS  THE  PROGRAM.’ 
PRINT  ♦ 

END  IF 

CONTINUE 

READdOUNIT,  12,END=100)  LINE 
FORMAT (A7 2 3 

IF(LINE.EQ. ’#’ )  STOP  '  *♦**  User  STOP  in  NAMELIST’ 

IF(LINE.EQ. ’?’ )  GOTO  1 
NEQ  =  INDEXfLINE, ’=’ ) 

IF(NEQ.EQ.O)  GOTO  100 
IF(LINE.EQ. ’QUIT’ )  THEN 
LINE=’QUIT=T’ 
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61  C 

62  C 

63  C 

64  C 

65 

66  13 

67 

68 

69 

70  C 

71 

72 

73 

74 

75 

76 

77 

78  15 

79  C 

80  C 

81 
82 

83  20 

84  C 

85 

86 

87 

88  25 

89  C 

90  C 

91 

92 

93  30 

94  C 

95 

96 

97  32 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108  35 

109  C 

110  C 

111 
112 

113  40 

114  C 

115 

116 

117 

118  C 

119  100 

120 


GOTO  13 

END  IF 
GOTO  100 
END  IF 

CONTINUE 

NAME=LINE(1:NEQ-1) 

NVAL=LINE(NEQ+1:72) 

DO  20  1=1, N1 

IF(NAME.EQ.R4NAMES(I))  THEN 
IFdNDEXINVAL.  ’  .  M.EQ.O)  THEN 

NPT  =  INDEXfNVAL, '  ') 

NVAL(NPT:NPT1  = 

END  IF 

READ(NVAL, 15)  R4VALUE(I) 

FORMAT! F15. 5) 

PRINT  *.  R4NAMES(I).'  =  ’.R4VALUE{I) 

PRINT  ♦ 

GO  TO  10 
END  IF 
CONTINUE 

DO  30  1=1, N2 

IF{NAME.EQ. I4NAMES(I>)  THEN 
READ(NVAL,25)  I4VALUE(I) 

FORMAT! 115) 

PRINT  *,  I4NAMES!!).’  =  M4VALUE(I) 

PRINT  *■ 

GO  TO  10 
END  IF 
CONTINUE 

DO  40  1=1, N3 

IF!NAME.EQ.L1NAMES!I))  THEN 
CONTINUE 

NDXT  =  INDEX! NVAL, 'T' ) 

NDXF  =  INDEX !N VAL, ’F' ) 

IF!NDXT*NDXF  .NE.  0  .OR.  !NDXT+NDXF) . EQ.  0)  THEN 
PRINT  *, NAME, ’IS  A  LOGICAL  VARIABLE.  ENTER  T  OR  F 
READ!6, 12)  NVAL 
GO  TO  32 
END  IF 

IF!NDXT.NE.O)  NVAL=’ . TRUE. ’ 

IF!NDXF.NE.O)  NVAL= FALSE . ’ 

READ! NVAL, 35)  LI VALUE! I) 

FORMAT !L15) 

PRINT  ♦,  LINAMES!!),’  =  ’,L1VALUE!I) 

PRINT  * 

GO  TO  10 
END  IF 
CONTINUE 

PRINT  ♦.’VARIABLE  NAME  ’.NAME,’  NOT  RECOGNIZED.’ 
PRINT  ♦. ’ INPUT  CONTINUES. . . ’ 

GO  TO  1 

CONTINUE 

PRINT  *,’User  input  complete.’ 
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121  PRINT  * 

122  RETURN 

123  END 


’DUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
'lUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


124 

:JUMBER  OF  WARNINGS  IN  COMPILATION  :  0 
'JUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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1 

2  C* 


3 

4 

5 

6 
7 
3 
9 

10 

11 

12  C****** 


C* 

C* 

C* 

C+ 

c+ 

c* 

c* 


H  A  R  R  Y  N  E  T 

A  RECOMFIGURABLE  30-NEURON  NETWORK  MODEL 
WHICH  LEARNS  BY  THE  DRIVE-REINFORCEMENT  LAW 


* 

+ 

* 


-.(c **+;*;***.'+,*■.*+• -ic^^c it, :(f  +  *  +  ;(c:(c 


SUBROUTINE  HARRYNETIT, KSROOT, KIROOT, KOROOT) 
INCLUDE  ’ \SYSPRO\COMNSH. INC’ 


13  C 


\SYSPRO\COMHSH .  INC  --  Abbreviated  labeled  common  arra?/£ 


for  use 


14 

C 

in  all  subroutines  except  EVOLVE. 

15 

C 

NEVER  CHANGE  ANYTHING  IN  THIS  FILE. 

16 

c 

Use  an 

INCLUDE  statement  to  use  these  common  arrays  in  any 

17 

C 

SYSPRO 

subroutine. 

18 

C 

19 

COMMON 

/STATSP/ 

STATEVC  1 ) 

20 

COMMON 

/KSNAME/ 

KSNAME! 2.  1) 

21 

COMMON 

/INPSP  / 

RINPUT!  1) 

22 

COMMON 

/KINAME/ 

KINAME! 2,  1) 

23 

COMMON 

/OUTPSP/ 

OUTPUT!  1) 

24 

COMMON 

/KONAME/ 

KONAME! 2,  1) 

25 

COMMON 

/OUTINT/ 

OUTINT! 2.  1) 

26 

COMMON 

/TIME  / 

TIME 

27 

CHARACTER*12  KSNAME, KINAME, KONAME, ISYSNM*6 

28 

COMMON 

/SIMVAR/ 

ENDTIM, MODE, DELTAT, TIMINC, NPRINT, AUDIT, RANDOM, 

29 

1 

NSYS.NXTSUB, ISYS!7, 110) , ISUB( 0:220) , ISYSNM! 110) , 

30 

1 

NPLOTS,  NSKIP,  KUR'/E!  5,  51 ) .  NPAGE,  RSMIN,  RSMAX,  RSEED 

31 

LOGICAL 

AUDIT. RANDOM 

32 

COMMON 

/  DTG  / 

ISEC, IMIN, IHR, IDAY, IMO,  lYR, 

33 

1 

JSEC. JMIN, JHR, JDAY, JMO, JYR, 

34 

1 

KSEC, KMIN, KHR, KDAY,  KMO, KYR 

35 

COMMON 

/TITLE  / 

ITITLE!40. 5) . IDATE, ITIME 

36 

CHARACTER 

ITITLE=K2.  IDATE*9,  ITIME*8 

37 

c 

38  C**** 

39 

40 

41 

42 

43  C 

44  0*^'  + 

45  C 

46  C* 

47  C 


**  END  OF  \SYSPRO\COMNSH. INC 
INCLUDE  ’ \BPNET\RUMDAT. INC’ 

COMMON  /RUMDAT/  GAMMA! 4 ) , PARMTHR( 5, 0 : 4 ) 

COMMON  /NETWORK/  NUMINPT,  INPUT! 50) , NUMNEUR, NEURN! 100) . 

NEDGE!2, 5100) . EDGEWT! 5100 ) , NFANIN! 4,  100) 


COMMON 


/WORKSP/  W!20,9) 


EXACTLY  AS  IN  MAIN  PROGRAM. 


C 

C 

C 


48 

49 

50 

51 

52 

53  C* 

54  C* 

55  C 


57 

58 


59  C* 

60 


VARIABLE  NAMING  SECTION 

**  -K  4:  He  ***■+■***♦*  * 

CHARACTER=t^l2  LSYSNM 

CHARACTER*12  T,SNAME!2,  1),  LINAME! 2 , 50 ) ,  LONAME!2,50) 
=  KSLEN  =  KILEN  =  KOLEN 

DIMENSION  KVNDX!  3,  0:80) 

THE  NUMBER  OF  SUBSYSTEMS.  N  =  NSUBS,  IS; 

V 

DATA  NSUBS  /  80/ 


DATA 
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62 

C 

THE  LENGTHS 

OF  THE  "SYSMODEL"  SYSTEM  VECTORS  (EXCLUDING 

63 

C 

SYSl  -  SYSN) 

ARE: 

64 

C 

65 

C* 

STATEV  : 

KSLEN 

vv  (  MUST  =  0  FOR  COMPOSITE  SYSTEM) 

66 

DATA 

KSLEN 

/ 

0/ 

67 

c* 

RINPUT  : 

KILEN 

=: 

vv 

63 

DATA 

KILEN 

/ 

50/ 

69 

c* 

OUTPUT  : 

KOLEN 

w 

70 

DATA 

KOLEN 

/ 

50/ 

71 

C 

72 

C 

THESE  VALUES 

ARE  REPORTED  TO  THE  CALLING  PROGRAM  BY  THE 

73 

74 

75 

76 

77 
73 

79 

80 
81 
82 
83 


84  C^t' 

85  C* 


86 

87 

88 

89 

90 

91 

92 

93 

94 


96 

97  C* 


98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 
109 

no 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 


"AUDIT"  SECTION  IF  AUDIT  =  .TRUE. 

IF  M>1  THEN  THIS  IS  A  COMPOSITE  SYSTEM  AND  IT  EVOLVES  THE 
STATEV  INDIRECTLY  BY  FIRST  EXECUTING  THE  CROSSTALK  FUNCTIONS 

(01 . ON)  TO  ADJUST  THE  INPUTS  TO  THE  SUBSYSTEMS  AND  THEN 

BY  CALLING  THE  SUBSYSTEM  MODELS  (SYSl . SYSN5. 

INTERMEDIATE  VALUES  (HOT  REQUIRED  TO  BE  KNOWN  UPON  ANY  ENTRY 
INTO  THIS  SYSMODEL  SUBROUTIHEi  SHOULD  BE  EQUIVALENCED  TO 
THE  WORKSPACE  VECTOR,  W  .  TO  SAVE  SPACE.  DO  THIS  NOW: 


EQUIVALENCE 


(  W(  1.1).  TEMPI  ) 
(ETC.  ) 


DATA 


NOTE : 


OUTPUT 

51. 


/ 


STATE  INPUT 
KVNDX/  1, 1. 1. 1.  51, 

237»0 

The  remaining  components  of  the  KVNDX  array  will  be  com¬ 
puted  in  the  AUDIT  SECTION  below,  on  the  assumption  that 
the  NEURON  subsystems  each  have  60  statevector  components, 
60  inputvector  components,  and  6  outputvector  components. 


SYSTEM  NAME: 
DATA  LSYSNM/ 


'KLOPF' 


THE  LENGTH  OF  THE  HISTORY  SEGMENT  FOR  EACH  SWAPSE  IS; 
DATA  LHIST/6/ 


STATEMENT  FUNCTION  SECTION 

NDS(J)  IS  THE  INDEX  OF  THE  J-TH  ENTRY  OF  THE  STATE  VECTOR 
OF  THIS  SYSTEM  (I.E.,  RELATIVE  TO  KSROOT).  AND  SIMILARLY  FOR 
NDKJ)  AND  NDO(J). 

NDS(J)  =  J  +  KSROOT 

NDKJ)  =  J  +  KIROOT 

NDO(J)  =  J  +  KOROCT 

NRS(I)  IS  THE  INDEX  OF  THE  ROOT  OF  THE  STATE  VECTOR  OF  THE 
I-TH  SUBSYSTEM  OF  THIS  SYSMODEL.  (ETC.  FOR  NRI.  NRO) 


NRS(I)  =  KSROOT  +  KVNDX(l.I)  -  1 
NRKI)  =  KIROOT  +  KVNDX(2,I>  -  1 
NRO(I)  =  KOROOT  +  KVNDX(3,I)  -  1 
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121  C 

122  C 

123  C 

124  C 

125  C 

126 

127 

128 

129 

130 

131  C 

132  C 

133  C 

134  C 

135  C 

136  C 

137  1000 

138 

139  C 

140 

141  C 

142 

143 

144 

145 

146 

147  1010 

148  C 

149  C 

150 

151 

152  1091 

153 

154  1093 

155 

156  1095 

157 

158  1096 

159  1020 

160  C 

161 
162 

163  C 

164  C 

165 

166  C 

167  C 

168  C 

169  C 

170  C 

171  2000 

172  C 

173  C 

174  C 

175 

176 

177 

178 

179 

180 


THE  FOLLOWING  STATEMENT  FUNCTIONS  SIMPLIFY  REFERENCES  TO  THE 
SYSTEM  VECTOR  ELEMENTS.  THEY  MAY  BE  USED  ONLY  ON  THE  RIGHT 
SIDE  OF  AN  ASSIGNMENT  STATEMENT.  ON  THE  LEFT  SIDE  OF  AN 
ASSIGNMENT  STATEMENT.  THE  FULL  REFERENCE  MUST  BE  USED. 

ST(I)  =  STATEV{NDS( I) ) 

Ria)  =  RINPUTfNDKI) ) 

OU(I)  =  OUTPUT(NDO(in 
STSUBCJ.I)  =  STATEV(NRS(J)+I) 

OUSUBfJ.I)  =  OUTPUTfNROIJ)+I) 

*  +  ***♦*:***♦♦:♦♦*  4  * 


AUDIT  SECTION 

4 +4+ 4444 *44*4 


CONTINUE 

IFI.HOT.  AUDIT)  GO  TO  2000 

IF(NUMNEUR.GT. NSUBS  .OR.  NUMNEUR. LT. 1 )  GOTO  5100 

NSUBS  =  NUMNEUR 
DO  1010  1=2, NSUBS 

KVNDXd.D  =  KVNDXd.l)  +  (I-D*  60 
KV^JDXd.I)  ^  KVNDXI2.1)  +  (I-l)t  63 
KVNDX(3.I)  =  KVNDX(3.1)  +  (I-l)*  7 
CONTINUE 

INITIALIZE  SYSTEM  VECTOR  LABELS 
DO  1020  J=l,50 
WRITE(LINAME(1. J) . 1091)  J 
FORMAK'NET  ’.12.’  INPUT’) 

WRITE(LINAME(2. J) . 1093) 

FORMAT! 'FIRING  RATE  ’) 

WRITEfLONAME! 1, J), 1095)  J 
FORMAT! ’NET  ’.12.’  OUTPT’ ) 

WRITE!LONAME!2, J) . 1096) 

FORMAT! ’SIGNAL  ’) 

CONTINUE 

CALL  SYSAUD! NSUBS. KSROOT. K I ROOT, KOROOT, KSLEN , KILEN . KOLEN , 
LSYSNM. LSNAME, LINAME, LONAME ) 

SKIP  THE  SUBSYSTEM  CROSSTALK  SECTION  DURING  AUDIT. 

GO  TO  3000 


SUBSYSTEM  CROSSTALK  SECTION 

4444444 4444444 4444444444444 

CONTINUE 

DISTRIBUTE  EXTERNAL  INPUTS  TO  THEIR  DESIGNATED  SYNAPSES. 

DO  2100  J  =  l.NUMINPT 
Y  =  RI!J) 
lEDGE  =  INPUT!J) 

N  =  NEDGE! 1. lEDGE) 

IF!N.EQ.0  .OR.  lEDGE.EQ.O)  GO  TO  2100 
DO  2050  K  =  lEDGE+l, lEDGE+N 
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181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 


NN  =  NEDGE(l.K) 

C  =  DESTINATION  NEURON  tt 

KSYN  =  NEDGE(2.K) 

C  =  DESTINATION  SYNAPSE  «  ON  DEST’N  NEURON 

IF (KSYN  .EQ.  0)  THEN 
OUTPUT(NRO(NN)+2)  =  Y 
ELSE 

C  SHIFT  PRIOR  INPUTS  TO  RIGHT 

JSYN  =  NRKNN)  +  ( KSYN-1 )  *(  LHIST+1 ) 

DO  2040  JL  =  LHIST, 1,-1 
JPL  =  JSYN  +  JL  +  1 
RINPUT(JPL)  =  RKJPL-l) 

2040  CONTINUE 

RINPUT( JSYN+1)  =  Y*EDGEWT(K^ 

END  IF 

2050  CONTINUE 
2100  CONTINUE 
C 

C  WARNING:  The  /NETWORK/  common  block  violates  the  system  simulation 
C  rules.  This  network  cannot  be  assembled  into  a  larger 

C  system.  Fold  it  into  /RINPUT/  before  trying  to  include 

C  BPNET  into  any  larger  SYSPRO  system. 

C 

C 

C  STATE  EVOLUTION  SECTION 

3000  CONTINUE 

C 

C 

DO  3010  J  =  l.NUMNEUR 

CALL  KLOPFON(T,NRS(J).NRI(J). NRO ( J ) ) 

3010  CONTINUE 

C 

IF(  AUDIT  )  RETURN 
C*  GO  TO  4000 

C 

C 

C  READOUT  SECTION 

C  + 

4000  CONTINUE 

C 

C  NOTE:  EACH  SUBSYSTEM  HAS  ITS  OWN  READOUTS.  THE  ONLY  READOUTS 

C  THAT  SHOULD  BE  INCLUDED  HERE  ARE  THOSE  THAT  USE  STATEV 

C  COMPONENTS  WHICH  ARE  NOT  ALL  MEMBERS  OF  A  SINGLE  SUB- 

C  SYSTEM  STATE  VECTOR. 

C 

C  DISTRIBUTE  EACH  NEURON’S  OUTPUTS  TO  THEIR  DESIGNATED  SYNAPSES 

C 

DO  4200  J  =  l.NUMNEUR 
lEDGE  =  NEURN(J) 

N  =  NEDGECl, lEDGE) 

IFdEDGE.EQ.O  .OR.  N.EQ.O)  GO  TO  4200 
DO  4150  K  =  lEDGE+l, lEDGE+N 
IF(NEDGE(1,K).NE.0)  THEN 
DO  4130  L=1,LHIST+1 
KL  =  NRI(NEDGE(1,K))+NEDGE(2.K)+L-1 
RINPUTCKL)  =  OUSUB(J. L)*EDGEWT(K) 

4130  CONTINUE 
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241 

ELSE 

242 

OUTPUT(NDO(NEDGE(2.K) ) 1  = 

OUSUB!J, 1) 

243 

END  IF 

244 

4150 

CONTINUE 

245 

4200 

CONTINUE 

246 

C 

247 

C 

248 

4999 

RETURN 

249 

C 

250 

C 

251 

C 

ERROR  RECOVERY  SECTION 

252 

C 

253 

5000 

CONTINUE 

254 

C 

255 

5100 

PRINT  5900,  NUMNEUR.NSUBS 

256 

STOP 

257 

5900 

FORMAT!'  SUBROUTINE  BPNET  - 

ERROR : 

YOU  HAVE 

’  .  14,  ’ 

258 

1 

'  ARRAY 

DIMENSIONS  ONLY 

ALLOW  ’ 

259 

C 

260 

END 

NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


NEURONS. ’ / 
.  14,  ’  .  ’  ,  /) 


NUMBER  OF  WARNINGS  IN  COMPILATION  :  0 
NUMBER  OF  ERRORS  IN  COMPILATION  ;  0 
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1  C* 

2  C* 

3  C* 

4  C* 

5  C* 

6  C* 

7  C* 

8  C* 

9  C* 

10 

11  c*** 

12 


♦  * 

*  KLOPFON  * 

*  * 

*  BASIC  PROCESSING  ELEMENT  MODEL  FOR  DRIVE-  ♦ 

*  REINFORCEMENT  NEURON  MODEL  FOR  BIOMASSCOMP  * 


**!f  !f*)f:*****#**if**)fr:f3f:*fr)f:#:f:)fsf*)f  if  ♦♦*♦**♦***;**  if 

JANUARY,  1988 

SUBROUTINE  KLOPFONd,  KSROOT,  KIROOT,  KOROOT) 

INCLUDE  ’\SYSPRO\COMNSH. INC’ 


13 

14 

15 

16 

17 

18 
19 


28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 


\SYSPRO\COMNSH. INC  —  Abbreviated  labeled  common  arrays,  for  use 

in  all  subroutines  except  EVOLVE. 

NEVER  CHANGE  ANYTHING  IN  THIS  FILE. 

Use  an  INCLUDE  statement  to  use  these  common  arrays  in  any 
SYSPRO  subroutine. 


20 

COMMON 

/STATSP/ 

STATEVf 

1) 

21 

COMMON 

/KSNAME/ 

KSNAME(2. 

1) 

22 

COMMON 

/INPSP  / 

RINPUT{ 

1) 

23 

COMMON 

/KINAME/ 

KINAME(2, 

1) 

24 

COMMON 

/OUTPSP/ 

OUTPUT { 

1) 

25 

COMMON 

/KONAME/ 

KONAME(2, 

1) 

26 

COMMON 

/OUTINT/ 

OUTINT(2. 

1) 

27 

COMMON 

/TIME  / 

TIME 

CHARACTER*12  KSNAME.  KINAME.  KONAME,  ISYSNMifS 

COMMON  /SIMVAR/  ENDTIM, MODE, DELTAT, TIMINC, NPRINT, AUDIT,  RANDOM, 

NSYS , NXTSUB , I SYS ( 7 , 1 1 0 ) , I SUB { 0 : 2  20 ) , I S YSNM (110), 
NPLOTS, NSKIP, KURVE( 5, 51 ) , NPAGE, RSMIN,  RSMAX,  RSEED 
AUDIT, RANDOM 

ISEC, IMIN, IHR,  IDAY,  IMO,  lYR, 

JSEC, JMIN, JHR, JDAY. JMO, JYR, 

KSEC, KMIN, KHR, KDAY,  KMO,  KYR 
ITITLE(40, 5), IDATE, ITIME 
ITITLE*2. IDATE*9,  ITIME*8 


LOGICAL 

COMMON 


/  DTG  / 


COMMON  /TITLE  / 
CHARACTER 


C**)f)f=fif  END  OF  \SYSPRO\COMNSH.  INC 
INCLUDE  ’\BPNET\RUMDAT. INC' 

COMMON  /RUMDAT/  GAMMA( 4 ) , PARMTHRf 5, 0 : 4 ) 

COMMON  /NETWORK/  NUMINPT, INPUT( 50) , NUMNEUR, NEURN( 100) , 

1  NEDGE(2,5100).EDGEWT(5100),NFANIN(4,  100) 

C 

C**if 


46  Cif 

47  C^ 

48 

49 

50  C 

51  C 

52  C 

53  C 

54  C 

55 

56 

57  C* 

58  C* 

59  C 

60 


COMMON  /WORKSP /  W ( 2  0 , 9 ) 

EXACTLY  AS  IN  MAIN  PROGRAM. 

DIMENSION  C(5) 

CHARACTER* 1  LETTER(0:5) 


VARIABLE  NAMING  SECTION 

ififififififififftifififififififififififififif 

CHARACTER* 12  LSYSNM 

CHARACTER*12  LSNAME{2. 60) ,  LINAME(  2, 63 ) ,  LONAMEd,  7) 
=  KSLEN  =  KILEH  =  KOLEN 

DIMENSION  KVNDX(  3.  1) 
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61  C* 

62  C 

63  C* 

64  C* 


EXTERNAL 


86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98  C 

99  C 
100 
101 
102 

103 

104 

105  C 

106  C 

107  C 


108  1000 
109 


C 

C 


110 
111 
112 

113 

114  1090 

115  C 

116 

117 

118 

119 

120 


THIS  NUMBER  MUST  =  MAX(l.NSUBS) 

DIFFEQ 

. . .  AND  ANY  OTHER  EXTERNAL  DECLARATIONS. 


65 

C 

66 

C 

THE  NUMBER  OF  SUBSYSTEMS,  N 

=  NSUBS,  IS: 

67 

c* 

V 

68 

DATA 

NSUBS  /  0/ 

69 

C 

70 

C 

THE  LENGTHS 

OF  THE  '"NEURON" 

SYSTEM  VECTORS 

(EXCLUDING 

71 

c 

SYS!  -  SYSN) 

ARE: 

72 

c 

73 

c* 

STATEV  : 

KSLEN  =  vv 

(TEN  SYNAPSES, 

6  HISTORICAL 

74 

DATA 

KSLEN  /  60/ 

75 

c* 

RINPUT  : 

KILEN  =  vv 

76 

DATA 

KILEN  /  63/ 

77 

c* 

OUTPUT  : 

:  KOLEN  =  w 

73 

DATA 

KOLEN  /  7/ 

79 

c 

80 

c 

THESE  VALUES 

ARE  REPORTED  TO 

THE  CALLING  PROGRAM  BY  THE 

81 

c 

"AUDIT"  SECTION  IF  AUDIT  = 

. TRUE. 

82 

c 

83 

DATA 

LETTER  /’a’ . 

’b’  ,  ’c’ ,  'd' ,  ’e’ 

,  'f  ’  / 

84 

DATA 

LHIST  /6/ 

85 

c 

STATEMENT  FUNCTION  SECTION 
^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

NDS(J)  IS  THE  INDEX  OF  THE  J-TH  ENTRY  OF  THE  STATE  VECTOR 
OF  THIS  SYSTEM  (I.E.,  RELATIVE  TO  KSROOT),  AND  SIMILARLY  FOR 
NDI(J)  AND  NDO(J). 

NDS(J)  =  J  +  KSROOT 

NDI(J)  =  J  +  KIROOT 

NDO(J)  =  J  +  KOROOT 


ST(I)  =  STATEV(NDS( I) ) 
RI(I)  =  RINPUT(NDI(I) ) 
OUd)  =  OUTPUT (NDO(  I) ) 


AUDIT  SECTION  (Calculations  needed  only  on  first  entry  can 
also  be  inserted  here.  ) 

CONTINUE 

IF(.N0T.  AUDIT)  GO  TO  2000 

INITIALIZE  THE  SYSTEM  VECTOR  LABELS 
NN  =  1+KSROOT/KSLEN 
WRITECLSYSNM,  1090)  NN 
FORMAK  ’PE  it’  .12) 

DO  1010  1=1,9 

IS  =  1  +  (I-1)*LHIST 
II  =  1  +  ( I-1)*(LHIST+1) 

DO  1009  J=0,LHIST 
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121 

122  1091 

123 

124  1092 

125 

126 

127  1093 

128  1094 

129 

130 

131 

132 

133  1095 

134 

135  1009 

136  1010 

137  C 

138  C 

139 

140 

141  C 

142  C 

143  C 

144 

145 

146  C 

147 

148 

149  1200 

150  C 

151 

152  C 

153 

154  C 

155 

156  C 

157 

158 

159  C 

160 

161  C 

162  C 

163  C 

164  2000 

165  C 

166  C 

167  C 

168  C 

169  C 

170  3000 

171  C 

172  C 

173  C 

174  C 

175  C 

176  C 

177 

178  C 

179  C 

180 


IF(J.LT.LHIST)  WRITE{LSNAME( 1. IS+J) . 1091)  I , LETTER( J ) . NM 
FORMAK’WT  Ml.Al.’  PE  ’.12) 

WRITElLINAMEf 1. II+J) , 1092)  I. LETTERS J) . NN 
FORMAK'INP  ’.Il.Al,’  PE  ’.12) 

IF(J. LT  LHIST)  WRITEfLSNAME(2. IS+J) , 1093) 

WRITEfLTMAME(2. II+J) . 1094) 

FORMAT (’D-R  SYNAPSE  ’) 

FORMAT ( ’ FREQUENCY  ’  ) 

IFd.EQ.  1)  THEN 

WRITE(LONAMEa.  I+J).  1095)  NN,  J+1 
END  IF 

FORMATf’PE  ’.12.’  OUT( ’ , 1 1 .  ’  )  ’  ) 

CONTINUE 

CONTINUE 


CALL  SYSAUDf NSUBS. KSROOT. KIROOT. KOROOT. KSLEN, KILEN, KOLEN, 

1  LSYSNM. LSNAME. LI NAME,  LONAME ) 

INITIALIZE  SOME  QUANTITIES  THAT  DON’T  CHANGE  FROM  ONE  NEURON 
TO  THE  NEXT: 

THETA  =  PARMTHR(l.O) 

=  FIRING  THRESHOLD  FOR  ALL  NEURONS 
DO  1200  I  =  1.5 
C(I)  =  PARMTHR(2.  I-l) 

CONTINUE 

=  LEARNING  RATE  CONSTANTS 
WMIN  =  PARMTHRO.O) 

=  LOWER  BOUND  FOR  ABSC SYNAPTIC  WEIGHTS) 

OrjMAX  =  PARMTHRM.O) 

=  UPPER  LIMIT  FOR  OUTPUT  LEVEL 
GAIN  =  PARMTHR(5.0) 

=  LEARNING  RATE  GAIN  FACTOR 


SKIP  THE  SUBSYSTEM  CROSSTALK  SECTION  IN  AUDIT  CYCLE 
GO  TO  3000 


CONTINUE 

SUBSYSTEM  CROSSTALK  SECTION 
<  NONE  FOR  PRIMITIVE  SYSTEM  > 


CONTINUE 

STATE  EVOLUTION  SECTION 

FIRST  TIME  THROUGH,  UNDO  SOME  OF  THE  (POSSIBLY  RANDOM) 
INITIALIZATION  OF  THE  STATE  VECTOR  WHICH  OCCURS  AFTER 
THE  AUDIT  SEQUENCE,  ERGO.  CAN’T  BE  DONE  ABOVE. 

IF(TIME. LT. TIMINC)  THEN 
END  IF 
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181 

182  C 

183  C 

184 

185 

186 

187  3100 

188 

189 

190 

191 

192 

193  C 

194  C 

195  C 

196  C 

197  C 

198 

199 

200 
201 

202  C 

203 

204 

205  C 

206  c 

207 

208  C 

209 

210  3200 

211 
212 

213  C 

214 

215 

216  3250 

217 

218  C 

219 

220 
221 
222 

223 

224 

225 

226  3300 

227 

228  C 

229  4000 

230  C 

231  C 

232  C 

233 

234  4100 

235 

236 

237 

238  C 

239  C 

240  C 


NN  =  1  +  KSROOT/KSLEN 

THIS  IS  THE  NUMBER  OF  THE  CURRENT  NEURON. 


Y  =  0. 

DO  3100  J  =  1.49,LHrST 

Y  =  Y  +  ST(J)*RI(J) 

CONTINUE 

Y  =  MINIOUMAX,  MAX(0. .  Y-THETA)) 
DY=  0U(1)-Y 


LEARNING  LAW 
♦  *+.**)(£>♦:*:*<**  + 

SYNAPTIC  LEARNING 

FOR  EACH  SYNAPSE  ( I )  ... 

DO  3300  1=  1.9 

KSY  =  LHIST*(I-1)+1 

KIN  =  (LHIST+1):<^(  I-l)  +  l 

INTEGRATE  OVER  THE  LAG  (J1  TO  OBTAIN  DELTA  W. 

DW=  0. 

DO  3200  J  =  1,LHIST-1 

(NOTE;  Should  go  to  LHIST,  but  that  v;ould  require  7  input 
lags,  i.e.,  one  more  than  there  are  synapse  lags.) 
DXIJ=  MAX(0. .  RI(KIN+J)-RI(KIN+J+1)  ) 

(this  implements  Klopf’s  refinement  on  p.  13) 

DW  =  DW  +  C(J)tABS(ST(KSY+J))*DXIJ 

CONTINUE 

DW  =  DY*DW»GAIN 

SHIFT  ALL  WEIGHTS  TOWARD  THE  PAST  (RIGHT  SHIFT) 

DO  3250  J  =  LHIST-l, 1,-1 

STATEV(NDS(KSY+J) )  =  ST(KSY+J-1) 

CONTINUE 

UPDATE  THE  CURRENT  VALUE  OF  THIS  SYNAPSE. 

WT  =  ST(KSY)  +  DW 
IF(ABS(WT).LT.  WMIN)  THEN 

STATEV(NDS(KSY) )  =  SIGN(WMIN, WT) 

ELSE 

STATEV(NDS(KSy) )  =  WT 
END  IF 

CONTINUE 


CONTINUE 
OUTPUT  SECTION 

DO  4100  J=LHIST+1, 2, -1 

OUTPUT(NDO(J) )  =  OU(J-l) 
OUTPUT(NDO( 1) )  =  Y 

RETURN 

ERROR  RECOVERY  SECTION 
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241  C 

242  5000 

243  C 

244  C* 

245  C* 

246 

247  C 

248 


CONTINUE 

INSERT  ERROR  RECOVERY  CODE  AND  MESSAGES  HERE. 

WRITE  TO  UNIT  5  (TERMINAL)  OR  UNIT  6  (LOGGING  FILE). 
GO  TO  3000 

END 


NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT;  0 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT;  0 


NUMBER  OF  WARNINGS  IN  COMPILATION  ;  0 
NUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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1  SUBROUTINE  HKDAT(T) 

2  C 

3  C  This  subroutine  provides  the  initial  simulation  parameters 

4  C  not  found  in  the  initialization  TRACE  file  and  provides  for 

5  C  the  computation  or  the  reading  of  the  time-varying  sensory 

6  C  inputs  to  the  neural  system. 

7  C 

8  INCLUDE  ’\SYSPR0\C0MNSH. INC’ 

9  c*>t'**** 


10 

C 

\SYSPRO\COMNSH.  INC  —  Abbreviated  labeled  common  arrays,  for  use 

11 

C 

in  all  subroutines  except  EVOLVE. 

12 

C 

NEVER  CHANGE  ANYTHING  IN  THIS  FILE. 

13 

C 

Use  an 

INCLUDE  statement  to  use  these  common  arrays  in  any 

14 

C 

SYSPRO 

subroutine. 

15 

C 

16 

COMMON 

/STATSP/ 

STATEVf  1) 

17 

COMMON 

/KSNAME/ 

KSNAME(2.  1) 

18 

COMMON 

/INPSP  / 

RINPUTI  1) 

19 

COMMON 

/KINAME/ 

KINAME(2.  1) 

20 

COMMON 

/OUTPSP/ 

OUTPUT  (•  1) 

21 

COMMON 

/KONAME/ 

KONAME! 2,  1) 

22 

COMMON 

/OUTINT/ 

0UTINT(2.  1) 

23 

COMMON 

/TIME  / 

TIME 

24 

CHARACTER=*fl2  KSNAME,  KINAME,  KONAME,  ISYSNMt6 

25 

COMMON 

/SIMVAR/ 

ENDTIM.MODE, DELTAT, TIMINC, NPRINT, AUDIT. RANDOM, 

26 

1 

NSYS.NXTSUB.  ISYS(7.  110),  ISUB(0:220),  ISYSNMdlO), 

27 

1 

NPLOTS, NSKIP, KURVE! 5. 51). NPAGE, RSMIN, RSMAX, RSEED 

28 

LOGICAL 

AUDIT, RANDOM 

29 

COMMON 

/  DTG  / 

ISEC, IMIN, IHR, IDAY, IMO, lYR, 

30 

1 

JSEC, JMIN, JHR, JDAY, JMO, JYR, 

31 

1 

KSEC, KMIN, KHR,  KDAY,  KMO.  KYR 

32 

COMMON 

/TITLE  / 

ITITLE(40. 5) . IDATE,  ITIME 

33 

CHARACTER 

ITITLE*2, IDATE*9, ITIME+8 

34 

C 

35  C*****»  END  OF  \SYSPR0\C0MNSH. INC 

36  INCLUDE  ’ \BPNET\RUMDAT. INC’ 

37  COMMON  /RUMDAT/  GAMMA(4) , PARMTHR(5, 0: 4) 

38  COMMON  /NETWORK/  NUMINPT,  INPUT(  50) ,  NUMIJEUR,  NEURN(  100) , 

39  1  NEDGE(2, 5100), EDGEWT(5100),NFANIN(4. 100) 

40  C 

41  C 

42  IF(T.GT. TIME)  GOTO  1000 

43  C 

44  D  PRINT  999,  ’SUBROUTINE  HKDAT:  READING  FILE  FORT2’ 

45  D  999  FORMAT (2 OX, A40) 

46  C 

47  C  READ  THE  SYSTEM  SIGMOID  PARAMETERS: 

48  READ(2,904) 

49  READ{2,903)  { (PARMTHR( J, I ) , 1=0. 4 ) . J=1 , 5 ) 

50  D  PRINT  903,  ( ( PARMTHR( J , I ) , 1=0, 4) , J=1 , 5 ) 

51  C 

52  C  READ  THE  NETWORK  GRAPH  STRUCTURE: 

53  READ(2,908) 

54  MAXINPUTS=50 

55  MAXNEURNS=80 

56  C  READ  THE  EXTERNAL  INPUT  DISTRIBUTION. . . 

57  C  ♦!t!*!t:*)(f***:»C!t:***************:tt*:)t******!(t 

58  NUMINPT  =  0 

59  lEDGE  =  0 

60  MIX  =  0 


RM /FORTRAN  Compiler  (V2.42) 
Source  File;  HKDAT.FOR 


Options;  /BLY 
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61  500 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85  C 

86  600 

87  C 

88  C 

89 

90 

91  610 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113  C 

114  700 

115  C 

116 

117 

118 

119  705 

120  710 


CONTINUE 

READ(2,907)  Ml. M2. M3. R4 

IF(M1.EQ.0  .AND.  M2.EQ.0  .A^JD.  M3.EQ.0)  GO  TO  600 

IFfR4. EQ. 0. )  R4=l. 

lEDGE  =  lEDGE  +  1 

IF(M1,EQ.M1X  .OR.  Ml.EQ.O)  THEN 

NEDGEf 1. INPUTfMlX) )  =  NEPGE{ 1 , INPUT(MIX) )  +  1 
NEDGEd.  lEDGE)  =  M2 
NEDGE(2. lEDGE)  =  M3 
EDGEWTdEDGE)  =  R4 
ELSE  IF(MI.GT.MIX)  THEN 

NUMINPT  =  MAX(NUMINPT.  MU 
INPUT(Ml)  =  lEDGE 
HEDGE ( 1. lEDGE)  =  1 
lEDGE  =  lEDGE  +  1 
HEDGE ( 1, lEDGE)  =  M2 
NEDGEf 2, lEDGE)  =  M3 
EDGEWT( lEDGE)  =  R4 
MIX  =  Ml 
ELSE 

PRINT  909 
STOP 

END  IF 
GO  TO  500 

CONTINUE 

READ  THE  INTERNAL  CONNECTION  GRAPH  STRUCTURE 

**)<£*:(c*4t)(r.;)i:*;t;>tt**=*<***>K**>t!****+;****  *********:*♦** 

READ (2. 904) 

MIX  =  0 

READ! 2, 907)  Ml.  M2. M3, R4 

IFIMl.EQ.O  .AND.  M2,E0..0  .AND.  M3.EQ.0)  GO  TO  700 

IF(R4.EQ.O. )  R4=l. 

lEDGE  =  lEDGE  +  1 

IFCMI.EQ.MIX  .OR.  Ml.EQ.O)  THEN 

HEDGE (l.NEURNI MIX))  =  NEDGE( 1, NEURN(MIX) )  +  1 
NEDGEI 1. lEDGE)  =  M2 
NEDGE(2. lEDGE)  =  M3 
EDGEWT( lEDGE)  =  R4 
ELSE  IFIMl.GT. MIX)  THEN 

NEURN(Ml)  =  lEDGE 
NEDGE( 1, lEDGE)  =  1 
lEDGE  =  lEDGE  +  I 
NEDGEf 1. I EDGE)  =  M2 
NEDGE(2. lEDGE)  =  M3 
EDGEWTf lEDGE)  =  R4 
MIX  =  Ml 
ELSE 

PRINT  910 
STOP 

END  IF 
GO  TO  610 

CONTINUE 

DETERMINE  THE  NUMBER  OF  NEURONS  IN  THE  NETWORK 
NUMNEUR  =  MIX 
DO  710  1=1, MIX 

DO  705  J=NEURN(I)+1,  NEURNf I ) +NEDGE( 1 , NEURNf I ) ) 
NUMNEUR  =  MAXfNUMNEUR, NEDGEf l.J) ) 

CONTINUE 
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121 

C 

122 

D 

123 

D 

124 

C 

125 

C 

126 

127 

128 

129 

130 

131 

132 

133 

C 

134 

135 

136 

D 

137 

138 

800 

139 

C 

140 

141 

805 

142 

143 

144 

145 

146 

147 

C 

148 

149 

C 

150 

1000 

151 

D 

152 

153 

154 

155 

156 

157 

158 

159 

160 

C 

161 

1100 

162 

163 

164 

C 

165 

901 

166 

902 

167 

903 

168 

904 

169 

905 

170 

906 

171 

907 

172 

908 

173 

909 

174 

175 

176 

910 

177 

178 

179 

911 

180 

1 

2 

1 

2 

1 


PRINT  *■,  'NUMBER  OF  INPUTS  =  ’  .  NUMINPT 
PRINT  *,  'NUMBER  OF  NEURONS^  ' . NUMNEUR 

CHECK  FOR  ERRORS  IN  THE  NETWORK  ARCHITECTURE 
IF(NUMINPT.GT.MAXINPUTS)  PRINT  9 1 1 , NUMINPT. MAXINPUTS 
IF(NUMNEUR.GT.MAXNEURNS)  PRINT  912, NUMNEUR. MAXNEURNS 
IF ( NUMNEUR. GT . MAXNEURNS . OR. NUMINPT . GT . MAXINPUTS )  STOP 
IF(NUMNEUR.LT.Ml)  THEN 

PRINT  *,  '  INPUT  ERROR.  HKDAT  COUNTED  TOO  FEW  NEURONS’ 

STOP 
END  IF 

READ(2, 904^ 

READ(2,902)  (GAMMA( 11 , 1=1 . 4) 

PRINT  902,  (GAMMAd),  1  =  1,4) 

lU  =  2 

CONTINUE 

SKIP  SIX  LINES  OF  INPUT  DATA  UNIT  (NEXT  READ  WILL  BE  ON  LINE  8) 
READ( lU, 906) 

CONTINUE 

READr lU, 901 . END=1 100 )  TT. N. G, Ml , R1 , M2 . R2 , M3 , R3 , M4, R4 
IF(TT.LT.O.  .AND.  IU.EQ.2)  THEN 
lU  =  Ml 
GO  TO  805 
END  IF 

RETURN 

CONTINUE 
PRINT  913,  TT 
IF(T.LT.TT)  RETURN 

IF(N.GT.O)  GAMMA(N)=G 
IF(Ml.GT.O)  RINPUT(M1)=R1 
IF(M2.GT.O)  RINPUT(M2)=R2 
IF(M3.GT.O)  RINPUT{M3)=R3 
IF(M4.GT.O)  RINPUT(M4)=R4 

READ( lU , 901 , END=1 100 )  TT, N, G , Ml , R1 . M2 . R2 , M3 , R3 . M4 . R4 
GO  TO  1000 

CONTINUE 
TT  =  l.OE+38 
GOTO  1000 

F0RMAT(F9. 3, IX. I 1 . IX, F8 . 3 , IX, 4( 13 , lX.FlO. 3, IX) ) 

FORMAT! 30X, 4F10. 8) 

FORMAT (20X, 5F10.  8) 

FORMAT!////) 

FORMAT! 20X, 5110) 

FORMAT!//////) 

FORMAT! 17. 18, 19, FIO. 4) 

FORMAT!////////) 

FORMAT! IX, 'SUBROUTINE  HKDAT  —  ERROR  READING  FILE  F0RT2 ’ / 

'  EXTERNAL  INPUTS  MUST  BE  LISTED  IN  ASCENDING  ORDER'/ 

'  EDIT  F0RT2  AND  RE-RUN  THE  PROGRAM’//) 

FORMAT! IX, 'SUBROUTINE  HKDAT  —  ERROR  READING  FILE  F0RT2 ’ / 

’  OUTPUT  NEURONS  MUST  BE  LISTED  IN  ASCENDING  ORDER’/ 

’  EDIT  F0RT2  AND  RE-RUN  THE  PROGRAM’//) 

FORMAT!’  SUBROUTINE  HKDAT  -  ERROR:  NUMBER  OF  INPUTS  =  ’,13./ 

’  MAXIMUM  »  INPUTS  =  ’,13.//) 


RM/FORTRAN  Compiler  (V2.42) 
Source  File:  HKDAT.FOR 


Options:  /BLY 
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181  912 

182 

183  913 

184 


FORMATf'  SUBROUTINE  HKDAT  ERROR:  NUMBER  OF  NEURONS  =  ’.13,/ 

’  MAXIMUM  tt  NEURONS  =  ’,13,//) 

FORMATf’  HKDAT;  READING  INPUT  DATA  FOR  TIME  =  ’.F9.3) 

END 


NUMBER  OF  WARNINGS  IN  PROGRAM  UNIT:  O 
NUMBER  OF  ERRORS  IN  PROGRAM  UNIT:  0 


NUMBER  OF  WARNINGS  IN  COMPILATION  :  0 
NUMBER  OF  ERRORS  IN  COMPILATION  :  0 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT:  BIOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  «:  1 

DATE:  1/11/88 

INPUT  CONDITIONS:  Two  random  variables  with  two  components 

each.  No  correlations  between  the  variables. 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

6000 

NX:  2 

NY: 

2 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN :  0 

X  ENTROPY: 

1.9574733 

Y  ENTROPY: 

1.9515425 

DMORPH: 

0. 0011546 

WHOLE  ENTROPY: 

3.9066575 

COMMENTS: 


FILE:  F0RT9;  ROUS;  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  — h 

X  ENTROPHV  X 
V  ENTROPHV  -  -B 


1.0  167.6  334.1  580.7  667.2  833.8  1000.3  1166.9  1333.4  1580.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  «: 
DATE: 


BIOMASSCOMP 

DltoRPH  Characterization 
David  G .  Boney 
2 

1/13/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  No  correlations  between  the  variables. 


VARIABLES: 


SEED:  .247E+13 

NX:  4 

LEXP: 

NY: 

3000 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN :  0 

X  ENTROPY: 

3.86821 

Y  ENTROPY: 

3.96083 

DMORPH: 

0.01346 

WHOLE  ENTROPY: 

7. 77527 

COMMENTS: 


FILE:  EXP2.DftT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -i — 

X  ENTROPHV  K 
V  ENTROPHV  -  -B 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 


RESEARCHER: 

EXPERIMENT  <t: 
DATE: 


David  G.  Boney 
3 

1/13/88 


INPUT  CONDITIONS:  Two  random  variables  with  six  components 

each.  No  correlations  between  the  variables. 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  6 

NY: 

6 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN:  0 

X  ENTROPY: 

5.87292 

Y  ENTROPY: 

5.89273 

DMORPH: 

COMMENTS : 

0. 30636 

WHOLE  ENTROPY: 

9.89462 

FILE:  EXP3.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL 


10.09  n 


vs. 


UH.  ENTROPHV 
X  ENTROPHV  K 

V  ENTROPHV  -  .e 
DftORPH- -  » 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  «: 
DATE: 


BIOMASSCOMP 

DMORPH  Characterization 
David  G.  Boney 
4 

1/13/88 


INPUT  CONDITIONS:  Two  random  vectors  with  eight  components 

each.  No  correlations  between  the  variables. 


VARIABLES: 


SEED:  .247E+13 

NX:  8 

LEXP: 

NY; 

3000 

8 

IC:  8 

A:  0.0 

B; 

1.0 

IFUN;  0 

X  ENTROPY: 

7.82463 

Y  ENTROPY: 

7.82070 

DMORPH: 

0. 53758 

WHOLE  ENTROPY: 

11.2658 

COMMENTS : 


FILE:  EXP4.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  — h 

X  ENTROPHV  X 
V  ENTROPHV  -  ^ 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


4 


i 

■1 


i 


t 


« 


« 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  »: 
DATE: 


BIOMASSCOMP 

DMORPH  Characterisation 
David  G .  Boney 
5 

1/15/88 


INPUT  CONDITIONS: 

Two 

random 

variables  with  four 

components 

each.  Intravariable  correlation. 

x(l)  =  x(2). 

VARIABLES: 

SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN :  1 

X  ENTROPY: 

3.8682065 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0.013411 

WHOLE  ENTROPY: 

7.7752690 

COMMENTS: 


FILE:  EXP5.DAT;  ROUS:  1  TO  1S00  :  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -+— 

X  ENTROPHV  K 
y  ENTROPHy  -  -e 
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PROJECT: 

EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT 

DATE: 

6 

1/15/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Two  intravariable  correlations,  x(l)  =  x(2>,  y(l)  =  y{2). 

VARIABLES: 

SEED:  .247E+13 
NX:  4 

A:  0. 0 

LEXP:  3000 

NY:  4  IC:  8 

B:  1.0  IFUN:  2 

X  ENTROPY: 

IM)RPH: 

3.8836992  Y  ENTROPY: 

0.0136949  WHOLE  ENTROPY: 

3.9413996 

7.7695165 

COMMENTS: 

FILE:  EXP6.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -f— 

X  ENTROPHV  X 
V  ENTROPHV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


« 


i 


I 


t 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT : 
EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  «: 
DATE: 


BIOMASSCOMP 

DMORPH  Characterization 
David  G.  Boney 
7 

1/15/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  y(l>. 

VARIABLES: 


SEED:  .247E+13 

NX:  4 

A:  0.0 

X  ENTROPY: 
DMORPH: 

COMMENTS : 


LEXP :  3000 

NY:  4 

B:  1.0 

3.9319389 
0. 2537242 


IC:  8 

IFUN :  3 

Y  ENTROPY:  3.9612269 

WHOLE  ENTROPY:  6.8684316 


FILE;  EXP7.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH,  ENTROPHV 

X  ENTROPHV  X 
y  ENTROPHV  -  .0 
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PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  #:  8 

DATE:  1/15/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l>  =  -y(l>. 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN:  4 

X  ENTROPY: 

3. 9063108 

Y  ENTROPY: 

3. 9612264 

DMOBPH: 

0.2536734 

WHOLE  ENTROPY: 

6.8430085 

COMMENTS: 


FILE:  EXP8.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -i— 

X  ENTROPHV  X 
V  ENTROPHV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT  tt: 
DATE: 

9 

1/15/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l>  =  .5  -  (yd)  -  .5). 

VARIABLES: 

SEED:  .247E+13 

NX:  4 

A:  0.0 

LEXP:  3000 

NY:  4 

B:  1.0 

IC:  8 

IFUN :  5 

X  ENTROPY: 

DMORPH: 

3.9063108 

0. 2536734 

Y  ENTROPY: 

WHOLE  ENTROPY: 

3.9612269 

6.8430085 

COMMENTS : 

FILE:  EXP9.DAT;  ROUS;  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV 

X  ENTROPHV  X 
V  ENTROPHV  -  -e 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Charac'teriza'tion 


RESEARCHER: 

EXPERIMENT  1»: 
DATE: 


David  G.  Boney 
10 

1/15/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation  xfl)  =  10  *  y{l>. 


VARIABLES: 


SEED: 

NX: 

.247E+13 

4 

LEXP: 

NY: 

3000 

4 

IC: 

8 

A: 

0.0 

B: 

1.0 

IFUN: 

6 

X  ENTROPY: 
DMORPH: 

3. 9319389 
0. 2537242 

Y  ENTROPY: 
WHOLE  ENTROPY: 

3. 9612269 
6.8684316 


COMMENTS: 


FILE:  EXP10.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -i— 

X  ENTROPHV  X 
V  ENTROPHV  - 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT:  BIOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 


RESEARCHER: 

EXPERIMENT  it: 
DATE: 


David  G.  Boney 
11 

1/15/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l>  =  .1  ♦  y( 1) . 

VARIABLES: 


SEED:  .247E+13 

NX:  4 

A:  0. 0 

LEXP: 

NY: 

B: 

3000 

4 

1.0 

1C:  8 

IFUN :  7 

X  ENTROPY: 
DMORPH: 

3.9319389 
0. 2537242 

Y  ENTROPY: 

WHOLE  ENTROPY: 

3.  9612269 
6.  8684316 

COMMENTS: 


FILE:  EXP11.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  ws.  UH.  ENTROPHV  -h— 

X  ENTROPHV  K 
V  ENTROPHV  -  ..0 


1.0  167.6  334.1  500.7  667. Z  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  «:  12 

DATE:  1/15/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l>  =  y(l)  +  y(2). 


VARIABLES: 


SEED: 

NX: 

A: 

.247E+13 

4 

0.0 

LEXP: 

NY: 

B: 

3000 

4 

1.0 

IC:  8 

IFUN:  8 

X  ENTROPY: 
DMORPH: 

3.8949640 
0. 1347252 

Y  ENTROPY: 

WHOLE  ENTROPY: 

3.9612269 
7.  3120666 

COMMENTS: 


FILE:  EXP12.DftT;  ROUS:  1  TO  1508  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -i — 

X  ENTROPHV  K 
V  ENTROPHV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1800.3  1166.9  1333.4  1500.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 


BIOMASSCOMP 

DMORPH  Characterization 


RESEARCHER: 

EXPERIMENT  #: 
DATE: 


David  G.  Boney 
13 

1/15/88 


INPUT  CONDITIONS:  Two  retndom  variables  with  four  components 


each. 

One  intercorrelation,  x(l) 

=  y(l>  *  y(2>. 

VARIABLES; 

SEED: 

NX: 

A; 

. 247E+13 

4 

0.0 

LEXP:  3000 

NY:  4 

B;  1.0 

IC:  8 

IFUN :  9 

X  ENTROPY; 
DMORPH: 

3.9062233 

0. 122994 

Y  ENTROPY: 

WHOLE  ENTROPY; 

3.9612269 

7.3707037 

COMMENTS: 

FILE;  EXP13.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -t— 

X  ENTROPHV  K 
V  ENTROPHV  -  -B 


I 


I 


i 


I 


!» 


e 


I 


I 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 


BIOMASSCOMP 

DMORPH  Characterization 
David  G.  Boney 


EXPERIMENT  Ht:  14 

DATE:  1/15/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Two  intervariable  correlations,  x(l)  =y{l)  +y(2>,  y(3>  = 
x(3)  +  x(4). 

VARIABLES: 


SEED:  .247E+13 

NX:  4 

A:  0. 0 

X  ENTROPY: 
DMORPH: 


LEXP:  3000 

NY:  4 

B:  1.0 

3.9319389 
0. 3654028 


IC:  8 

IFUN:  10 

Y  ENTROPY: 
WHOLE  ENTROPY: 


3. 9216778 
6. 367x360 


COMMENTS: 


FILE:  EXP14.DAT;  ROUS;  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV 

X  ENTROPW  X 
V  ENTROPHV  - 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 


RESEARCHER: 

EXPERIMENT  »: 
DATE: 


David  G.  Boney 
15 

1/15/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  oomponents 

each.  Two  intervariable  correlations,  x(l)  =  y{l>  *  y(2),  y(3)  = 
x(3)  ♦  x(4). 

VARIABLES; 


SEED:  .247E+13  LEXP; 

NX:  4  NY: 

A;  0.0  B: 


3000 

4 

1.0 


IC;  8 

IFUN:  11 


X  ENTROPY: 
DMORPH: 


3. 9062233 
. 2310253 


Y  ENTROPY:  3.9316173 

WHOLE  ENTROPY:  6.8979411 


COMMENTS : 


FILE;  EXP15.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPHV  -i— 

X  ENTROPHV  K 
V  ENTROPHV  -  .0 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  Dt*3RPH  Characterization 

RESEARCHER;  David  G.  Boney 

EXPERIMENT  «:  16 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Two  intervariable  correlations,  x(l)  =  y(l),  x(2)  =  y(2). 
VARIABLES: 

SEED:  .247E+13  LEXP:  3000 

NX:  4  NY;  4  IC:  8 

A:  0.0  B:  1.0  IFUN;  12 

X  ENTROPY:  3.9312069  Y  ENTROPY:  3.9612269 

DMORPH:  0.4973724  WHOLE  ENTROPY;  5.8836598 

COMMENTS; 


FILE:  EXP16.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV 


K  ENTROPV  X 
V  ENTROPV  -  -a 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  17 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Three  intervariable  correlation,  x(l)  =  y(l),  x(2>  =  y(2), 
x(3)  =  y(3). 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN:  13 

X  ENTROPY: 

3.8895742 

Y  ENTROPY: 

3. 9612269 

DMORPH: 

0. 7411187 

WHOLE  ENTROPY; 

4. 8575912 

COMMENTS : 


FILE:  EXP17.DAT;  ROUS;  1  TO  1500  ;  PLOT  OF.  TRIAL  vs.  UH.  ENTROPV 


X  ENTROPV  X 
V  ENTROPV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  #: 

DATE: 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Four  intervariable  correlations,  x(l>  =  y(l),  x(2)  =  y(2) 
x(3)  =  y(3),  xf4)  =  yM). 

VARIABLES: 


BIOMASSCOMP 

DMORPH  Characterization 
David  G.  Boney 
18 

1/18/88 


SEED: 

NX:  4 

A:  0.0 


247E+13  LEXP: 
NY: 


X  ENTROPY; 
DMORPH: 

COMMENTS; 


B: 

3. 9612269 
0.9807996 


3000 

4 

1.0 


IC: 

IFUN: 


8 

14 


Y  ENTROPY: 
WHOLE  ENTROPY: 


3. 9612269 
3.9612269 


FILE:  EXP18.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL 


3.95 

3.52 

3.08 

2.64 

2.20 

1.76 

1.32 

0.88 

0.44 

0.00 


vs.  UH.  ENTROPV  -h — 

X  ENTROPV  . X 

V  ENTROPV  -  -0 
^ - l—a - ?< - a - K - 8 - 1 - Jf-B - k  DMORPH - H  . j*.. 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

trial 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT:  BIOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 


RESEARCHER: 


David  G .  Boney 


EXPERIMENT  «:  19 

DATE:  1/18/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

Ooe  intervariable  correlation,  x(l>  =  .1  *  y(i)  +  9  * 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

1C:  8 

A:  0.0 

B: 

1.0 

IFUN:  15 

X  ENTROPY: 

3.8948736 

Y  ENTROPY: 

3.9612269 

DMORPH: 

COMMENTS: 

0. 2156875 

WHOLE  ENTROPY: 

6.9849877 

FILE:  EXP19.DAT:  ROUS;  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV 

X  ENTROPV  . .K- 

V  ENTROPV  -  -a 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT:  BIOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  #:  20 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  .2  *y(l)  +  .8  * 
y(2). 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN:  16 

X  ENTROPY: 

3.8984380 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0. 1874058 

WHOLE  ENTROPY: 

7.  1027756 

COMMENTS: 


FILE:  EXP20.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  ws.  UH.  ENTROPV  -t— 

K  ENTROPV  X 
V  ENTROPV  -  ^ 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT  «: 
DATE: 

21 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four 
each.  One  intervariable  correlation,  x(l)  =  .3  *  yd) 
y(2). 

VARIABLES: 

components 
+  .  7  * 

SEED:  .247E+13 

NX:  4 

A:  0.0 

LEXP:  3000 

NY:  4 

B:  1.0 

IC:  8 

IFUN:  17 

X  ENTROPY: 

DMORPH: 

3.8973515 

0. 1651944 

Y  ENTROPY: 

WHOLE  ENTROPY; 

3. 9612269 
7. 1913958 

COMMENTS ; 

FILE:  EXP21.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -i — 

X  ENTROPV  X 
V  ENTROPV  - 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G .  Boney 

EXPERIMENT  #: 
DATE: 

22 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  .4  *  yd)  +  .6  * 
y(2). 

VARIABLES: 

SEED:  .247E+13 

NX:  4 

A:  0.0 

LEXP:  3000 

NY:  4 

B:  1.0 

IC:  8 

IFUN:  18 

X  ENTROPY: 
DMORPH: 

3.8918359 

0.  1431046 

Y  ENTROPY:  3.9612269 

WHOLE  ENTROPY:  7.2750959 

COMMENTS: 

FILE:  EXP22.DAT;  ROUS: 

1  TO  1500  ;  PLOT  OF  TRIAL 

vs.  UH.  ENTROPV  -i— 

X  ENTROPV  K 
V  ENTROPV  -  -B 


! 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  # :  23 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  .5  *y^l)  +  .6  * 
y(2). 

VARIABLES: 


SEED:  .247E+13 

NX:  4 

LEXP: 

NY: 

3000 

4 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN:  19 

X  ENTROPY: 

3.8949640 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0. 1347252 

WHOLE  ENTROPY: 

7.3120666 

COMMENTS: 


FILE'  EXP23.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  — i 

X  ENTROPV  K 
V  ENTROPV  -  -B 


1.0  167.6  334.1  500. 7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT  K: 
DATE: 

24 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  y(l>  +  y(2)  +  y{3). 

VARIABLES: 

SEED:  .247E+13 

NX:  4 

A:  0.0 

LEXP:  3000 

NY:  4 

B:  1.0 

TP  •  ft 

IFUN:  20 

X  ENTROPY: 

DMORPH: 

3.9074244 

0. 1411841 

Y  ENTROPY: 

WHOLE  ENTROPY: 

3.9612269 
7. 2984409 

COMMENTS: 

FILE:  EXF24.DAT;  ROUS; 

1  TO  1500  ;  PLOT  OF  TRIAL 

vs.  UH.  ENTROPV  - 

j _ 

X  ENTROPV  X 
V  ENTROPV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1080.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  «:  25 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x(l)  =  (  y(l)  +  y( 2)  +  y( 3) 
)  /  3. 

VARIABLES: 


SEED;  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0.  0 

B; 

1.0 

IFUN:  21 

X  ENTROPY: 

3. 9074244 

Y  ENTROPY: 

3. 9612269 

DMORPH: 

0. 1411841 

WHOLE  ENTROPY: 

7. 2984409 

COMMENTS; 


FILE:  exp25.dat:  ROUS;  1  TO  1S00  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV 

X  ENTROPV  V 
y  ENTROPV  - 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT : 
EXPERIMENT: 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT  H: 
DATE: 

26 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  One  intervariable  correlation,  x^l)  =  y{ 1)  +  y(2)  +  y(3)  + 
y(4)  . 

VARIABLES: 

SEED:  .247E+13 
NX:  4 

A:  0.0 

LEXP;  3000 

NY:  4 

B:  1.0 

IC:  8 

IFUN:  22 

X  ENTROPY: 

DMORPH: 

3.8960431 

0. 1365768 

Y  ENTROPY: 
WHOLE  ENTROPY: 

3. 9612269 
7.3056674 

COMMENTS: 

FILE;  EXP26.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -+— 

X  ENTROPV  K 
V  ENTROPV  -  -a 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT : 
EXPERIMENT : 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT  K: 
DATE: 

27 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four 

each.  One  intervariable  correlation,  x(l)  =  Ty(l>  + 
y(4))  /  4. 

VARIABLES: 

components 
yI2>  +  y(3)  + 

SEED:  .247E+13 

NX:  4 

A:  0.  0 

LEXP:  3000 

NY:  4 

B:  1.0 

IC:  8 

IFUN:  23 

X  ENTROPY: 

DMORPH; 

3.8960431 

0. 1365768 

Y  ENTROPY: 

WHOLE  ENTROPY; 

3.9612269 

7.3056674 

COMMENTS : 

FILE:  EXP27.DAT;  ROUS;  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -i — 

X  ENTROPV  K 
y  ENTROPV  -  -a 
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APPENDIX  C 

DMOEPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 


BIOMASSCOMP 

DMORPH  Characterization 
David  G .  Boney 


EXPERIMENT  #:  29 

DATE:  1/18/88 


INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Two  intervariable  correlations,  x(l)  =  y(l)  +  y{2),  x(2)  = 
yl2)  +  y{3). 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

A;  0.0 

B: 

1.0 

X  ENTROPY: 

3.8991964 

DMORPH: 

0. 2496906 

IC:  8 

IFUN:  25 

Y  ENTROPY:  3.9612269 

WHOLE  ENTROPY;  6.8519797 


COMMENTS: 


FILE:  EXP29.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV 

X  ENTROPV  K 
V  ENTROPV  - 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


Page  C. 28 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT : 
EXPERIMENT : 

BIOMASSCOMP 

DMORPH  Characterization 

RESEARCHER: 

David  G.  Boney 

EXPERIMENT 

DATE; 

30 

1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four 

each.  Three  intervariable  correlations,  xfl)  =  y(l)  + 
=  y{2)  +  y(3),  x(3)  =  y(3)  +  y(4). 

VARIABLES: 

components 
y(2),  x^2) 

SEED;  .247E+13 
NX:  4 

A;  0. 0 

LEXP;  3000 

NY:  4 

B :  1.0 

IC:  8 

IFUN:  26 

X  ENTROPY; 

DMORPH: 

3.8101032 

0.  3473946 

Y  ENTROPY: 

WHOLE  ENTROPY: 

3. 9612269 
6. 3682823 

COMMENTS: 

FILE:  EXP30.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPy 

X  ENTROPV  X 
V  ENTROPV  -  -0 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT ;  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  it:  31 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Four  intervariable  correlations,  x(l)  =y(l)  +y{2),  x(2)  = 
y(2)  +  y(3).  x(3)  =  y(3)  +  y(4).  x(4)  =  y(4)  +  yd). 

VARIABLES: 


SEED:  .247E+13 

LEXP: 

3000 

NX:  4 

NY: 

4 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN:  27 

X  ENTROPY: 

3.8055966 

Y  ENTROPY: 

3.9612269 

ZMIRPH: 

0. 4399206 

WHOLE  ENTROPY: 

5. 9900842 

COMMENTS: 


FILE:  EXP31.DAT;  ROUS:  1  TO  1S00  ;  PLOT  OF  TRIAL  ws.  UH.  ENTROPV  -i — 

X  ENTROPV  K 
V  ENTROPV  -  -0 . 

5.95 
5.29 

4.63 

3.97 

3.31 

2.64 

1.98 

1.32 
0.66 
0.00 

1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 


APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  D^K)RPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  tt:  33 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Two  intervariable  correlations.  xfl>  =  y(l)  +  y(2)  +  y(3), 
x(2)  =  y(2>  +  y(3)  +  y(4). 

VARIABLES: 


SEED:  .247E+13 

LEXP; 

3000 

NX;  4 

NY; 

4 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN:  29 

X  ENTROPY: 

3.9173765 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0. 2494148 

WHOLE  ENTROPY; 

6.8712735 

COMMENTS: 


FILE;  EXP33.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV 


X  ENTROPV  X 
y  ENTROPV  -  -B 


1.0  167.6  334.1  500.7  667.2  833.8  1008.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT:  BIOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  It:  34 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each. 


VARIABLES : 


SEED:  .247E+13 

NX:  4 

LEXP: 

NY: 

3000 

4 

IC:  8 

A:  0.0 

B: 

1.0 

IFUN:  30 

X  ENTROPY: 

3.9007394 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0.3511075 

WHOLE  ENTROPY: 

6. 4439230 

COMMENTS: 


FILE:  EXP34.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -+— 

X  ENTROPV  X 
V  ENTROPy  -  -a 


1.0  167.6  334.1  500.7  667.2  833.8  1008.3  1166.9  1333.4  1500.0 

TRIAL 
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APPENDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT ;  B lOMASSCOMP 

EXPERIMENT:  DMORPH  Characterization 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  <»:  35 

DATE:  1/18/88 

INPUT  CONDITIONS:  Two  random  variables  with  four  components 

each.  Four  intervariable  correlations  x^l)  =  y(l)  +  y(2)  +  y(3), 
x(2)  =  y(2>  +  y(3)  +  yM).  x(3)  =  y(3)  +  yf4)  +  ya).  x(4)  =  y(4) 
+  y{l)  +  y(2). 

VARIABLES: 


SEED:  .247E+13 

NX:  4 

LEXP: 

NY: 

3000 

4 

IC:  8 

A:  0. 0 

B: 

1.0 

IFUN:  31 

X  ENTROPY: 

3.7893291 

Y  ENTROPY: 

3.9612269 

DMORPH: 

0. 4477465 

WHOLE  ENTROPY: 

5.9422097 

COMMENTS : 


FILE:  EXP35.DAT;  ROUS:  1  TO  1500  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -h — 

K  ENTROPV  X 
V  ENTROPV  -  ■€ 


1.0  167.6  334.1  500.7  667.2  833.8  1000.3  1166.9  1333.4  1500.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  -  BACK  PROPAGATION 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  1»:  1 

DATE:  2/1/88 

INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-3-4  back  propagation 
network  that  is  suppose  to  pass  through  its  inputs.  The  four 
components  are  the  following  sin  functions:  x( 1)  =  .5  +  .5  ♦ 
sin(t)  ,  x(2)  =.54.5*  sin(t  -  1>,  x{3>  =.54  .25  *  sin(t), 
x(4)  =  .5  4  .25  *  sin(t-l).  t  is  the  simulation  time.  The  y 
variable  is  a  vector  of  four  components  that  is  the  output  of  the 
network.  The  run  was  done  with  learning  off  and  the  output  was 
sampled  at  the  bound ry  of  a  second. 

VARIABLES: 

LEXP:  1500 

NX:  4  NY:  4  IC:  90 

X  ENTROPY:  3.1543148  Y  ENTROPY:  3.5262272 

DMORPH:  0.3564391  WHOLE  ENTROPY:  5.0859146 

COMMENTS : 
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PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  -  BACK  PROPAGATION 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  #:  2 

DATE:  2/1/88 

INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-3-4  back  propagation 
network  that  is  suppose  to  pass  through  its  inputs.  The  four 
components  are  the  following  sin  functions:  x(l)  =  .5  +  .5  ♦ 
sin(t)  ,  x(2)  =  .5  +  .5  *  sin(t  -  1>,  x(3>  =  .5  +  .25  ♦  sin(t), 
x{4)  =  .5  +  . 25  *  sin(t-l).  t  is  the  simulation  time.  The  y 
variable  is  a  vector  of  four  components  that  is  the  output  of  the 
network.  The  run  was  done  with  learning  off  after  having  learned 
for  1500  seconds.  Ssunpling  was  done  at  the  second  boundries. 
VARIABLES: 

LEXP:  1500 

NX:  4  NY:  4  IC:  90 

X  ENTROPY:  3.1543148  Y  ENTROPY:  3.1566122 

DMORPH:  .3363257  WHOLE  ENTROPY:  4.6819711 

COMMENTS: 


Page  C. 35 


APPEHDIX  C 

DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  -  BACK  PROPAGATION 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  it :  3 

DATE:  2/1/88 

INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-3-4  back  propagation 
network  that  is  suppose  to  pass  through  its  inputs.  The  four 
components  are  the  following  sin  functions:  x(l)  =  .5  +  .5  * 
sin(t)  ,  x(2)  =  .5  +  .5  *  sin(t  -  1),  xC3>  =  .5  +  .25  ♦  sin(t), 
x(4)  =  . 5  +  . 25  *  sin{t-l).  t  is  the  simulation  time.  The  y 
variable  is  a  vector  of  four  components  that  is  the  output  of  the 
network.  This  run  was  done  with  learning  off  and  sampled  once  a 
second  at  the  half  second  boundry. 

VARIABLES: 

LEXP:  1500 

NX;  4  NY:  4  IC:  90 

X  ENTROPY:  3.1540511  Y  ENTROPY:  3.3335972 

DMORPH:  0.4456712  WHOLE  ENTROPY:  4.4079671 

COMMENTS: 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT :  B lOMASSCOMP 

EXPERIMENT:  DMORPH  -  BACK  PROPAGATION 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  4 

DATE:  2/1/88 

INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-3-4  back  propagation 
network  that  is  suppose  to  pass  through  its  inputs.  The  four 
components  are  the  following  sin  functions:  x(l>  =  .5  +  .5  * 
sin(t>  ,  x(2)  =  .5  +  .5  ♦  sin(t  -  1),  x(3)  =  .5  +  .25  *  sin(t), 
x(4>  =  .5  +  .25  ♦  sin(t-l).  t  is  the  simulation  time.  The  y 
variable  is  a  vector  of  four  components  that  is  the  output  of  the 
network.  This  run  was  done  with  learning  off  after  having  run  for 
1500  seconds  with  learning  on.  the  sampling  was  done  once  a 
second  on  the  half  second  intervals. 

VARIABLES: 


LEXP:  1500 

NX:  4 

NY:  4 

IC;  90 

X  ENTROPY: 

3. 1540511 

Y  ENTROPY; 

3. 1570110 

DMORPH: 

0. 5151951 

WHOLE  ENTROPY: 

3.8159778 

COMMENTS; 
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PIMJJECT ;  B  lOMASSCOMP 

EXPERIMENT:  DMORPH  -  KLOPF 

RESEARCHER:  David  G.  Boney 

EXPERIMENT  #:  1 

DATE:  1/29/88 

INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-2-4  Klopf  network. 

The  four  components  are  the  following  sin  functions:  x(l>  =  .5  + 
.5  *  sin(t)  .  xf2)  =  .5  +  .5  *  sin{t  -  1),  x(3)  =  .5  +  .25  * 
sin(t),  x(4)  =  .  5  +  .  25  ♦  sin(t-l).  t  is  the  simulation  time.  The 
y  variable  is  a  vector  of  four  components  that  is  the  output  of 
the  network.  The  inputs  and  outputs  where  in  three  second 
intervals.  This  run  was  done  with  learning  off 
VARIABLES: 


LEXP:  1475 

NX:  4 

NY:  4 

IC:  90 

X  ENTROPY: 

3. 1604311 

Y  ENTROPY: 

3.6817343 

DMORPH: 

0.4795595 

WHOLE  ENTROPY: 

4. 7713003 

COMMENTS: 

FILE;  HK4BASE.JU9;  ROUS:  1  TO  1475  ;  PLOT  OF  TRIAL  vs.  UH.  ENTROPV  -i— 


X  ENTROPV  K 
V  ENTROPV  -  -0 


1.0  164.8  328.6  492.3  656.1  819.9  983.7  1147.4  1311.2  1475.0 

TRIAL 
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DMORPH  EXPERIMENTS  AND  GRAPHS 


PROJECT: 

EXPERIMENT: 

RESEARCHER: 

EXPERIMENT  it: 
DATE: 


BIOMASSCOMP 
DMORPH  -  KLOPF 

David  G.  Boney 

2 

1/29/88 


INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-2-4  Klopf  network. 

The  four  components  are  the  following  sin  functions:  x(l)  =  .5  + 
.5  *  sin(t>  ,  x(2)  =  .5  +  .5  *  sin(t  -  1),  x(31  =  .5  +  .25  * 
sin(t),  xM)  =  .5  +  .25  ♦  sin(t-l).  t  is  the  simulation  time.  The 
y  variable  is  a  vector  of  four  components  that  is  the  output  of 
the  network.  The  inputs  and  outputs  where  in  three  second 
intervals.  This  run  was  done  with  learning  off  after  having  run 
for  1500  seconds  with  learning  on. 

VARIABLES: 


LEXP:  1475 

NX:  4 

X  ENTROPY: 
DMORPH: 

COMMENTS : 


NY:  4 

3. 1604311 
. 4120096 


IC:  90 

Y  ENTROPY:  3.1796041 

WHOLE  ENTROPY:  4.3539858 


FILE;  HK4AFTR.JU9;  ROUS;  1  TO  1475  ;  PLOT  OF  TRIAL  vs.  UH.  EHTROPV  -i— 

X  ENTROPy  v 
y  ENTROPy  -  ^ 


1.0  164.8  328.6  492.3  656,1  819.9  983.7  1147.4  1311.2  1475.0 

TRIAL 
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PROJECT: 
EXPERIMENT : 

RESEARCHER; 

EXPERIMENT  #; 
DATE: 


BIOMASSCOMP 
DMORPH  -  KLOPF 

David  G.  Boney 

3 

1/29/88 


INPUT  CONDITIONS:  The  x  variable  is  a  vector  with  four 

components.  This  vector  is  the  input  to  a  4-2-4  Klopf  network. 

The  four  components  are  the  following  sin  functions:  x(l)  =  .5  + 
.5  ♦  sin(t)  ,  x{2)  =  .5  +  .5  *  sin(t  -  1).  x(3)  =  .5  +  .25  * 
sin(t),  x(4)  =  . 5  +  . 25  *  sin(t-l).  t  is  the  simulation  time.  The 
y  variable  is  a  vector  of  four  components  that  is  the  output  of 
the  network.  The  inputs  and  outputs  where  in  three  second 

after  changing  two  of  the  neuron 

coefficients,  running  the  network  for  1500  seconds  with  learning 
on,  and  then  running  for  1500  seconds  with  learning  off  The 
sampling  was  done  with  leaning  off. 

VARIABLES: 


LEXP:  1475 

NX;  4  NY:  4 


IC:  90 


X  ENTROPY;  3.1604311 

DMORPH:  .4120096 


Y  ENTROPY;  3.1176403 

WHOLE  ENTROPY:  4.3872814 


FILE:  HK4AFTR2.JU9;  ROUS:  1  TO  1475  ;  PLOT  OF  TRIAL  vs.UH.  ENTROPV 

X  ENTROPV  K 
V  ENTROPV  -  ^ 


1.0  164.8  328.6  492.3  656.1  819.9  983.7  1147.4  1311,2  1475.0 

TRIAL 
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ABSTRACT 


Realtime  signal  processing  of  mul ti el ectrodc  probes  of  living 
neural  networks  is  limited  both  by  the  speed  and  flexibility  of 
the  host  computing  equipment  and  by  the  efficiency  and  flexibility 
of  the  signal  processing  algorithm.  In  this  report,  we  describe 
the  structural  and  functional  design  of  a  multichannel  signal 
processing  algorithm  which  has  the  ability  to  dynamically  include 
or  exclude  processing  steps  and  subroutines  in  order  to  maximize 
th€->  utilization  of  available  hardware.  That  is,  the  algorithm 
will  perform  all  the  processes  that  it  can  perform  on  realtime 
data  in  a  racetrack  buffer  without  either  overtaking  the  incoming 
data  or  falling  behind  and  being  lapped  thereby.  Remaining 
prctcesses  are  performed  on  intermediate  stared  data  in  an  off-line 
(non-realtime)  mode. 


I 
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1.  INTRODUCTION 

1.1  Review  o-f  the  Problem 

The  primary  problem  iddressed  by  this  research  is  the  -fact 
that  it  is  now  possible  far  most  neurophysiology  laboratories  to 
collect  more  parallel  channels  o-f  data  -from  living  neuron  networks 
than  they  can  af-ford  to  either  save  or  to  process  in  real  time. 

For  e>!ample,  Dr.  Gross's  MMEP  apparatus  is  now  collecting  34 
parallel  channels  of  data  at  20K  samples  per  second  each.  Simply 
transferring  that  data  at  2  bytes  per  sample  to  tape  or  disk  for 
subsequent  non-realtime  processing  would  fill  up  a  50MB  volume  in 
a  little  over  half  a  minute. 

Even  though  a  half  minute's  worth  of  data  is  enough  to 
provide  useful  information  on  the  network  structure,  the  MMEP 
apparatus  has  the  potential  to  not  only  listen  to  the  culture 
network,  but  to  talk  back  to  it  as  will.  Without  realtime 
analysis  of  the  current  patterns  in  the  network,  this  valuable 
potential  cannot  be  exploited.  That  is,  while  it  would  be 
possible  to  "shout"  at  the  network  at  random  intervals  and  then 
analyze  the  reactions  offline  later,  the  truly  earthshaking 
experiments  that  could  be  performed  require  the  detection  of 
developing  patterns  of  signals  in  the  network  (in  real  time)  and 
the  selective  feedback  of  signals  to  interrupt  or  respond  to  those 
patterns. 

With  the  capability  for  realtime  interaction,  it  would  be 
possible  for  the  first  time  to  test  certain  mathematical  models  of 
neural  network  behavior,  such  as  the  synaptic  plasticity  models 
which  are  used  to  formulate  explicit  mechanisms  for  the  Hebbian 
learning  laws  (cf.,  Grossberg,  or  Hestenes  [21).  Tn  particular, 
Grossberg's  "outstar  learning  theorem"  could  be  tested  by 
repeatedly  injecting  a  pattern  of  signals  to  coincide  with  the 
occurrence  of  a  pattern  in  a  naturally  occurring  sequence.  Later, 
if  the  artificial  stimulus  evokes  the  same  response  as  the  natural 
pattern,  but  in  the  absence  of  the  natural  pattern,  and  in  the 
absence  of  any  prior  ability  to  evoke  that  response,  then  the 
required  synaptic  plasticity  will  have  been  demonstrated.  This 
could  have  enormous  consequences  for  science,  with  implications 
not  only  for  neurophysiology,  but.  for  psychology  and  computer 
science  as  well. 

To  do  the  necessary  realtime  processing  on  which  these 
experiments  are  predicated  —  indeed,  to  even  analyze  the 
culture's  network  behavior  offline  in  nonrealtime  —  requires  the 
development  of  mathematical  tools,  computational  algorithms,  and 
hardware  configurations  that  may  not  now  exist.  Since  the 
hardware  selection  is  limited  more  by  economic  considerations  than 
by  technological  capabilities,  we  are  driven  to  try  to  coax  the 
resulting  algorithms  to  be?  "hardware  friendly".  That  is,  we  want 
the  algorithm  to  be  able  to  adapt  to  the  host  without  excessive 
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reprogrammi  nq .  The  meaning  o-f  this  statement  will  become  more 
clear  in  section  1.3  (Objectives);  however,  we  shall  now  review 
the  existing  techniques  for  analysis  of  multichannel  neural 
network  data. 

1.2  Selection  of  Net hods 

The  tasks  which  must  be  performed  by  the  computer  in  order  to 
analyze  the  network  data  can  be  summarized  in  the  following 
sequence; 


(1)  1  Detect  spikes  in  the  signals  which  1 

1  are  sensed  in  each  electrode.  ! 


V 


(2)  !  Classify  the  detected  spikes  ! 

I  according  to  their  source  neurons.  ! 


(3)  I  Compress  the  data  streams  from 

1  each  source  neuron  into  a  minimal 
!  stream  still  containing  enough 
1  information  to  un-compress  and 
!  recover  the  original  data. 


V 


(4)  !  Process  the  compressed  data  from 

{  all  source  neurons  to  identify  the 
!  communication  structure  of  the 
I  total  network 


Tlte  steps  in  this  sequence  become  more  difficult  as  one  proceeds, 
until  at  step  (4)  one  finds  almost  nothing  beyond  some  rather 
straightforward  histograms  being  attempted  in  the  current  litera¬ 
ture.  Since  the  histograms  are  informative,  we  shall  provide  for 
their  computation,  but  we  shall  also  attempt  a  more  networl;- 
thE'oretic  description  of  the  system. 
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1.2.1  Spike  Detection 

Spike  detection  is  clearly  the  easiest  step  to  perform.  One 
"merely"  sets  a  threshold  at  a  level  which  excludes  the  noise,  and 
whenever  the  <absolute  value  of  the)  digitized  voltage  level 
exceeds  the  threshold,  one  declares  that  a  spike  has  been 
detected.  We  put  "merely"  in  quotes  because  the  selection  of  the 
threshold  level  in  a  multichannel  process  will  be  done  by  the 
computer  according  to  some  algorithm,  and  this  algorithm  is 
considerably  less  transparent  than  the  thresholding  instruction. 

Establishing  the  threshold  —  an  easy  task  for  a  person 
looking  at  an  oscilloscope  trace  —  is  essentially  a  problem  in 
"constant  false-alarm  rate"  (CFAR)  methods.  Digital  CFAR  thresh¬ 
olding  is  discussed  by  Rohling  in  C5D.  In  the  absence  of  spiking 
signals,  i.e.,  when  only  noise  is  being  sensed,  the  threshold  can 
be  determined  by  computing  a  histogram  of  the  digitized  voltages 
and  finding  a  level  within  which  enough  of  the  samples  lie,  so 
that  only  once  in  a  specified  number  of  seconds  does  a  noise 
sample  lie  outside  the  threshold. 

Unf or tunatel y  we  usually  have  to  take  the  noise  as  it  comes: 
with  some  signal  added  to  it.  Therefore  in  order  to  set  the  CFAR 
threshold  we  have  to  mask  out  the  spikes  so  that  they  are  excluded 
from  the  histogram.  Of  course,  we  can't  use  thresholding  to  find 
the  spikes,  because  it  is  the  threshold  we  are  trying  to  find!  So 
we  have  to  use  some  other  a-priori  knowledge  to  find  and  mask  the 
spikes.  If  the  computational  application  of  this  knowledge  were 
quicker  thc\n  thresholding,  then  we  would  naturally  use  it  for  the 
realtime  spike  detection;  but  it  isn't,  so  we  apply  it  prior  to 
realtime  processing  in  a  so-called  "learning"  mode. 

There  are  a  number  of  different  masking  techniques,  depending 
on  the  application,  but  one  that  comes  to  mind  exploits  the 
continuity  (smoothness)  of  the  spike  trace  as  opposed  to  the 
roughness  of  the  noise.  Thus  the  algorithm  will  examine  a  number 
of  consecutive  samples  to  see  if  their  magnitudes  are  all 
unusually  large  (compared  to  an  unmasked  histogram)  and  all  on  the 
same  side  of  the  origin.  If  so,  an  appropriate  amount  of  data 
around  the  suspected  spike  is  masked  off,  i.e.,  removed  from  the 
histogram.  This  process  continues  to  be  performed  on  the  learning 
data  set  until  no  further  spikes  can  be  identified,  whereupon  the 
remaining  histogram  represents  an  approximation  of  the  probability 
density  function  of  the  noise  alone.  The  threshold  can  then  be 
selected  at  that  value  for  which  the  number,  n,  of  samples  whose 
magnitudes  exceed  the  threshold,  divided  by  the  length,  T,  of  the 
learning  samiple  in  seconds  is  most  nearly  equal  to  the  desired 
false  alarm  rate. 

One  could  detect  spikes  by  algorithms  that  are  more 
complicated  than  simple  thresholding,  such  as  by  matched 
filtering,  maximum  likelihood  estimation,  and  many  others  (each 
requiring  its  own  "learning  mode"  and  each  also  requiring  its  own 
thresholding  operation),  but  these  methods  can  be  reserved  for 
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some  post-detection  processing  in  the  remaining  steps  of  the 
sequence. 

1.2.2  Spike  Cl assi f i cati on 

The  next  task  is  to  determine  whether  the  signal  on  a  given 
electrode  is  the  superposition  of  spike  trains  from  more  than  one 
neuron  source,  and  if  so,  to  separate  the  signal  into  its 
constituent  parts.  A  survey  of  computer  methods  for  this 
separation  task  was  given  by  Schmidt  C6]  in  1984.  Other 
techniques  can  also  be  found  <cf,  Okada  and  Maruyama  C43  ).  Our 
approach  in  this  section  is  to  select  a  sequence  of  methods  based 
on  a  progression  from  the  computationally  simple  to  the  more 
difficult.  We  choose  this  approach  because  the  structural  design 
of  our  algorithm  (Chapter  2)  calls  for  the  use  of  simple  methods 
on  channels  where  simple  methods  suffice,  and  harder  methods  where 
they  are  required,  thus  using  available  processing  time  most 
efficiently.  (Hardware  configurations  employing  a  dedicated 
signal  processor  in  each  channel  will  not  need  to  avail  themselves 
of  this  choi ce. ) 

As  with  the  detection  process,  spike  classification  is  done 
in  two  parts.  During  a  learning  mode  each  channel  is  evaluated  to 
determine  the  values  of  some  variables  which  will  be  used  in  the 
real  time  mode  and  which  are  expected  to  change  very  slowly,  if  at 
all.  For  spike  separation,  these  variables  will  partition  the 
channels  into  subsets,  each  of  which  can  be  de-i nter 1 eaved  by  a 
different  class  of  algorithm. 

The  first  of  these  variables  will  identify  the  number  of 
sources  being  received  on  the  channel.  If  only  one  source  is 
being  received,  then  step  (2)  of  the  sequence  is  trivial.  The 
second  variable  will  be  a  vector  whose  components  identify  the 
amplitudes  at  which  distinct  sources  are  found  and  the  number  of 
distinct  sources  that  are  found  at  each  amplitude.  If  no 
amplitude  bin  has  more  than  one  source  contributing  to  it,  then 
the  classification  can  be  performed  by  amplitude  discrimination. 
But  if  any  bin  has  more  than  one  source  then  some  form  of  waveform 
discrimination  will  have  to  be  used  to  separate  the  sources. 

The  easiest  way  to  determine  that  there  are  more  than  one 
source  neuron  represented  in  the  signal  is  to  perform  a  one¬ 
dimensional  cluster  analysis  on  the  peak  voltage  in  each  detected 
spike.  (A  general  description  of  cluster  analysis  algorithms  is 
given  in  Appendix  A.)  If  more  than  one  cluster  is  found,  then 
there  are  more  than  one  source.  The  fact  that  the  converse  of 
that  statement  is  false  necessitates  the  use  of  more  complicated 
algorithms  for  the  classification  process,  but  since  these  more 
complicated  mthods  work  best  when  the  peaks  are  presented  in 
constant-amplitude  clusters,  we  may  as  well  do  the  easiest  job 
first. 
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To  do  the  amplitude  (peak)  cluster  analysis,  one  must  first 
detect  and  measure  the  height  of  the  greatest  local  extremum 
within  a  single  pulse-width  of  the  threshold  crossing  (detected 
in  step  1).  this  can  generally  be  done  without  any  multiplies  or 
divides.  (See  Appendix  B,  or  reference  C4D).  Using  a  cluster 
diameter  that  is  large  enough  to  allow  for  known  amplitude 
variations  from  single  sources,  one  then  sorts  the  peak  values 
into  a  histogram  (see  Tou  &  Gonzalez  163  >  to  find  the  clusters. 

To  determine  whether  the  peaks  assigned  to  a  given  cluster 
break  down  further  into  different  wave  shapes  one  can  then  do  an 
n-di mensi onal  vector  cluster  analysis,  where  n  is  the  number  of 
samples  after  a  threshold  crossing  needed  to  cover  all  spike 
waveforms.  The  n-dimensional  clusters  are  found  analogously  to 
the  way  the  amplitude  clusters  are  found,  but  the  distance  measure 
is  a  little  more  involved  and  the  cluster  diameter  is  more 
difficult  to  establish  (See  Appendix  A). 

Once  the  waveform  clusters  are  identified,  it  might  be 
possible  to  do  the  realtime  assignment  of  spikes  to  the 
appropriate  cluster  with  an  algorithm  that  is  simpler  than  a 
template  comparison  or  a  matched  filter.  But  even  with  a 
"hardware  friendly"  algorithm,  it  is  fair  to  assume  that  there  is 
a  vector  or  array  processor  lurking  in  a  wait-state  nearby, 
eagerly  contemplating  its  next  victim.  Therefore,  we  shall  prefer 
to  simply  vectorize  the  distance  between  the  realtime  peak  and  the 
centers  of  the  clusters  to  assign  it  to  its  source. 


1.2.3  DATA  COMPRESSION 

The  first  step  in  data  compression  is  almost  taken  care  of 
in  tasks  (1)  and  (2)  simply  by  detecting  and  classifying  the 
spikes  in  each  channel.  By  delivering  a  report  of  the  TIME  when 
the  spike  was  detected  over  the  threshold  the  SOURCE  which  emitted 
it,  and  (perhaps)  the  measured  spike  amplitude,  one  has  compressed 
the  20  or  so  sample  values  representing  the  spike,  and  the  80  or 
so  preceding  sample  values  representing  noise  alone  into  only  one 
or  two  values  from  which  a  replica  of  the  spike  (sans  noise)  can 
be  reproduced. 

One  essential  item  in  the  description  of  the  SOURCE  (though 
it  may  actually  be  irrelevant  as  far  as  the  neuronal  "message"  is 
concerned)  is  the  wave  shape  of  the  spike,  which  was  discovered  in 
the  classification  step.  Since  this  shape  is  expected  to  change 
only  very  slowly,  if  at  all,  during  the  experiment,  it  can  be 
identified  (via  pointers  or  links  to  a  template  library)  with  the 
array  into  which  the  TIME  values  are  reported.  Thus,  if  one  wants 
to  resurrect  a  replica  of  the  raw  data  which  was  recorded  from  a 
particular  source  neuron,  one  can  retrieve  the  sequence  of  times- 
of -arrival  of  spikes  from  the  array  associated  with  that  source, 
and  at  those  times  construct,  a  pulse  with  the  shape  in  the 
template  library  that  is  linked  to  the  source.  If  amplitudes  are 
considered  important,  the  (normalized)  templates  can  be  scaled  by 
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the  amplitude  that  was  stored  alongside  the  time  oi  arrival. 

Data  compression  is  possible  whenever  there  are  replicated 
patterns  in  the  data  that  can  be  parametrized  and  reduced  to  a 
symbol,  -followed  by  a  list  o-f  values  o-f  the  parameters.  Thus,  in 
the  previous  paragraph,  raw  voltage  samples  -from  an  electrode 
which  is  sensing  two  sources  is  compressed  into  the  -following 
symbols  and  parameter  lists: 

SOURCE 1 ((T11,A11),<T12,A12),  ...) 

S0IJRCE2(  (T21,A21  )  ,  <T22,A22)  ,  .  .  .  )  , 

where  Tij  is  the  time  o-f  arrival  of  the  j-th  pulse  from  the  i-th 
source,  and  Ai j  is  its  amplitude.  The  name,  SOURCEl,  is  taken  to 
be  equivalent  to  the  unchanging  characteristics  of  the  source. 
Similarly,  mathematicians  use  symbols  like  "SIN",  and  "LOG"  to 
compress  the  descriptions  of  families  of  functions,  and  by 
supplying  parameters  like  the  frequency  of  the  sinusoid,  and  the 
base  of  the  logarithm,  they  can  then  resurrect  a  graphical 
representation  of  the  function. 

"Bursting"  i s  an  observable  feature  of  signals  drawn  from 
certain  kinds  of  neurons.  Even  though  the  literature  shows  little 
agreement  on  a  definition  of  what  might  constitute  a  burst,  we  can 
sidestep  that  issue  for  the  purpose  of  data  compression.  For  our 
purposes,  it  is  sufficient  to  establish  a  library  of  certain 
patterns  occurring  in  the  compressed  spike  reports,  and  when  those 
patterns  are  detected,  reduce  them  to  a  more  compact 
representation.  We  suggest  the  following  definition; 

A  "burst"  is  either  <1)  three  or  more  consecutive  spikes,  whose 
amplitude  sequence  lies  within  a  martingale  envelope  rooted  on  the 
first  spike,  and  whose  interval  sequence  lies  within  a  martingale 
envelope  rooted  on  the  first  interval;  or  (2)  any  single  spike 
which  fails  to  fit  in  a  sequence  of  the  previous  category.  The 
envelope  to  be  used  on  the  amplitude  sequence  should  be  of  the 
form. 


a  <b— a )  e>:p  ( -kt )  +/-  vt  , 

where  the  values  of  a,  b  and  k  can  be  determined  from  the 
first  three  amplitudes  in  the  candidate  sequence  (since  vt  is 
nearly  zero  at  the  start),  and  v  is  the  variance  of  the 
amplitude  process.  Similarly,  the  envelope  to  be  used  on  the 
interval  sequence  should  be  of  the  form, 

b’  -  mt  -*-/-  v't  , 

where  b'  and  m  can  be  determined  from  the  first  two  intervals, 
and  v'  is  the  variance  of  the  interval  process. 

Thus,  in  a  burst,  one  expects  the  amplitudes  to  fall  off 
along  some  declining  exponential  starting  at  the  first  pulse,  plus 
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or  minus  a  margin  that  gets  a  little  wider  as  the  burst  proceeds; 
and  one  expects  the  interpulse  spacing  to  increase  approximately 
linearly  with  time,  again  with  a  margin  that  increases  down  the 
line.  The  parameters  a,  b,  b’,  k,  and  m  are  measured  in  a  learn 
mode,  and  taken  to  be  characteristic  oP  the  source;  while  v  and 
v'  are  taken  to  be  rejection  thresholds  outside  oP  which  the 
candidate  burst  sequence  terminates. 


1.2.4  NETWORK  ANALYSIS 

The  state  oP  the  art  in  network  communication  analysis  oP 
living  neuron  networks  is  still  rather  primitive,  which  is  to  be 
expected  due  to  the  quite  recent  emergence  oP  the  technology  Por 
simultaneous  sensing  oP  numerous  points  within  isolated  networks. 
The  principal  tool  Par  the  analysis  is  the  pairwise  correlation  oP 
Peatures  oP  the  spike  trains  by  way  oP  the  cross-correlation 
histogram.  These  methods  are  described  in  Chapter  10  oP  MacGregor 
and  Lewis  C33,  and  in  several  other  papers  <e.g.,  C13). 

We  Peel  that  although  these  correlograms  are  usePul  tools  Por 
sparsely  connected  networks  such  as  might  be  Pound  in  aplesia  or 
in  sensory  ganglia,  we  are  not  likely  to  Pind  statistically 
signiPicant  correlations  appearing  in  the  more  densely  connected 
networks.  The  reasoning  here  is  that  networks  such  as  Pound  in 
the  mammalian  cortex  are  structured  Por  the  ePPicient  sorting  oP 
coordinated  patterns  oP  input  signals,  rather  than  Por  serving  as 
in-line  ampliPiers  oP  single  inputs.  Consequently,  it  is  only 
when  a  synaptic  input  is  a  part  oP  a  coordinated  pattern  oP  inputs 
that  it  will  participate  in  the  generation  oP  spiking  or  bursting 
at  the  output  oP  the  aPPerent  neuron. 

In  order  that  these  experiments  with  the  MMEP  apparatus 
should  Pul  Pi  11  their  potential,  we  Peel  that  they  should  result 
not  in  the  publication  oP  a  report  that  is  Pull  oP  histograms  and 
other  statistical  humdrum,  but  rather  that  they  should  be  used  to 
conPirm  or  eliminate  specipic  quantiPiable  hypotheses  regarding 
the  possible  mechanisms  oP  learning,  recall,  synaptic  plasticity, 
memory  storage,  and  the  like.  To  do  this  requires  the  use  oP 
models  which  link  the  hypotheses  to  certain  parameters  which  are 
susceptible  to  measurement  with  the  apparatus. 

One  such  model  is  provided  in  the  system  oP  coupled  nonlinear 
ordinary  diPPerential  equations  known  as  Grossberg's  Field 
Equations  (nicely  presented  by  Hestenes  in  C23).  These  equations 
describe  the  incremental  ePPect  on  the  ionic  potential  energy  oP 
certain  pulse  patterns  arriving  at  the  synapses.  The  details  are 
important,  but  they  can  be  summarized  by  pointing  out  that  Por  the 
excitatory  synapses,  the  ePPect  oP  the  incoming  signal  in  driving 
the  neuron  toward  its  Pirinq  threshold  is  proportional  to  the 
recent  history  oP  the  pulse  repetition  Prequency  (F'RF)  on  that 
synapse.  The  constant  oP  proportionality  is  a  characteristic  oP 
the  synapse  that  can  be  modi  Pied  by  the  coincidence  oP  neuronic 
Piring  and  input  to  the  synapse.  It  is  called  the  synaptic 
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coupling  coe-f -f  icient .  The  amount  of  history  that  is  relevant  is 
controlled  by  the  rate  at  which  the  neuron  will  dissipate  its 
energy  without  firing. 

According  to  this  model,  one  would  expect  to  find 
correlations  between  the  onset  of  bursts  from  a  given  neuron  and 
the  PRF  history  on  its  synapses.  This  suggests  that  converting 
the  spike  data  to  the  corresponding  PRF  history,  and  "stacking" 
these  histories  whenever  a  burst  occurs  at  the  output  of  a 
particular  neuron,  might  accumulate  strong  peaks  on  signals  that 
are  connected  to  excitatory  synapses.  On  the  other  hand,  if  the 
stacks  are  triggered  at  the  onset  of  blank  intervals  in  the 
neuron’s  output,  then  the  presence  of  strong  peaks  would  indicate 
an  inhibitory  influence.  In  contrast  to  the  cross-correl at i on 
(interval)  histogram,  which  attaches  significance  to  the  influence 
of  a  single  spike  at  the  input  to  a  subsequent  single  spike  at  the 
output,  this  technique  attaches  significance  to  the  recent  history 
of  ionic  current-pumpi ng  to  the  onset  of  bursting.  The  assumption 
here  is  that  the  subsequent  spikes  in  a  burst  are  a  "ringing" 
effect  due  to  the  close  coupling  of  each  neuron  to  itself,  rather 
than  a  direct  effect  of  the  input  signals.  Therefore,  if  they 
were  used  as  reference  points  for  stacking  the  input  PRF  signals, 
they  would  only  contribute  noise  and  computational  burden. 


1.3  SOFTWARE  DEVELOPMENT  OBJECTIVES 

The  structural  design  of  the  signal  processing  software  that 
is  presented  later  in  section  2  is  guided  by  the  following 
considerations.  First  is  the  need  to  make  efficient  use  of  the 
processing  hardware  in  the  MASSCOMP  5700  so  that  the  least  amount 
of  available  data  from  the  MMEP  apparatus  is  lost.  Preliminary 
timing  calculations  have  shown  that  so  long  as  the  preprocessing 
tasks  (tasks  1  and  2  from  page  3>  must  be  handled  by  the  MASSCOMP 
it  will  not  be  possible  to  do  any  burst  detection  or  network 
analysis  in  real  time,  and  probably  only  the  top  6  to  10  channels 
in  the  priority  list  can  be  reduced  to  spike  data.  The  remaining 
analysis  tasks  will  have  to  be  done  off-line  in  non-realtime  using 
previously  saved  spike  data. 

However,  the  next  consideration  is  that  hardware  development 
is  under  way  for  the  offloading  of  the  preprocessi ng  from  the 
MASSCOMP  to  an  array  of  TMS  32020  processors.  Therefore,  the 
processing  software  needs  to  be  flexible  in  its  scheduling  of 
processing  tasks,  according  to  whether  its  data  is  coming  directly 
from  the  MMEP  in  raw  form,  indirectly  from  the  MMEP  through  the 
32020  boards  as  spike  data  (but  still  in  real  time),  or  directly 
from  disk  or  tape  storage  as  spike  data  or  burst  data  on  demand. 

The  fact  that  the  software  is  being  developed  in  response  to 
experimental  necessities  requires  careful  attention  to  structured 
programming  and  modular  design.  Our  initial  expectations  of  the 
signal  processing  routines  that  will  be  effectual  for  the  desired 
analyses  will  have  to  be  modified  with  experience.  The  ability  to 
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hang  new  and  different  subroutines  into  the  scheduler  without 
incurring  timing  or  synchronization  problems  that  propagate 
unchecked  through  the  rest  of  the  processes  must  be  built  in  from 
the  start. 

Because  of  the  limitations  on  the  size  and  the  intensity  of  j 

the  programming  effort,  it  is  recognized  that  it  will  not  be 
feasible  to  attempt  to  implement  a  large  and  complex  software 
system.  On  the  other  hand,  it  is  highly  likely  that  a  low  level 
of  programming  effort  will  be  applied  to  the  project  over  a  number 
of  years.  It  is  wise,  then,  to  take  a  lesson  from  this 
investigator’s  past:  In  1980  I  was  assigned  to  restructure  a 

system  analysis  program  that  was  written  to  predict  the  | 

performance  of  a  large  solar  photovoltaic  energy  system.  That 

program  was  begun  small  and  grew  as  the  system  developed.  It  was 

in  the  form  of  a  single  FORTRAN  main  program  without  a  single 

subroutine  call  aside  from  intrinsic  functions!  It  was  several 

thousand  lines  of  incomprehensible  rat’s  nest.  I  am  not 

suggesting  that  anyone  on  this  project  would  be  quite  that  crude.  . 

Rather,  I  am  emphasizing  that  it  is  all  right  to  design  a  modular,  ' 

structured,  and  comprehensive  signal  processing  program  that  is 

perhaps  overwhelming  in  its  scope,  but  that  i s  at  least  not  likely 

to  have  to  be  razed  several  years  down  the  road.  With  that 

thought  in  mind,  we  proceed  to  the  structural  design  of  the 

algorithm. 

I 

2.0  STRUCTURAL  DESIGN 

The  structural  design  of  the  processing  software  is  specified 
by  the  HIPO  ("Heirarchy,  Input-Process-Output")  charts  which  are 
included  in  Appendix  C  to  this  report.  The  following  paragraphs 

are  intended  to  elaborate  on  those  charts  to  assist  the  ! 

programming  team  in  their  implementation. 

2.1  TOP  LEVEL  HIPO  DESCRIPTION 

The  top  level  HIPO  chart  contains  five  primary  processes. 

The  first  process  controls  the  initialization  of  the  program  and  j 

its  parameters  to  reflect  the  experimental  configuration  and  the 

proper  disposition  of  output  data  (filenames,  display  devices, 

etc,).  The  user  also  specifies  his  processing  priorities  to 

override  defaults  that  will  be  accepted  by  the  scheduler. 

The  second  process  accepts  data  from  the  selected  source 
devices  or  data  files,  and  performs  "learn  mode"  operations  on  it  i 

in  non-real  time.  These  operations  provide  the  preliminary 
pattern  recognition  functions  to  establish  thresholds,  identify 
clusters  and  cluster  centers,  and  tailor  a  processing  sequence  to 
each  channel  for  the  benefit  of  the  master  scheduler. 

The  third  process  accepts  data  from  the  designated  source  and 
applies  the  appropriate  data  compression  subroutines  in  accordance 
with  the  proce^ssing  requirements  and  timing  limitations.  This 
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process  does  not  analyze  the  compressed  data  (except  insofar  as 
some  form  of  analysis  is  inherent  in  the  compression),  but  instead 
provides  various  levels  of  compressed  data  to  facilitate  the 
analysis  routines  in  the  fourth  process.  (The  functions  to  be 
performed  are  described  in  paragraph  1,2  and  in  Appendices  A  and 
B.  ) 

The  fourth  process  applies  certain  statistical  and  analytical 
functions  in  accordance  with  user  specifications  and  timing 
limitations.  This  process  obtains  data  from  the  various  levels  of 
compression  for  graphical  representation,  listings,  etc.;  it 
computes  histograms,  correl ograms,  stacks;  and  it  provides 
transfer  of  compressed  data  and  processing  results  to  appropriate 
output  files/devices.  This  process  provides  the  primary  user 
interaction  with  the  ongoing  experiment. 

The  fifth  process  terminates  the  experiment.  It  is 
responsible  for  purging  buffers,  closing  files,  appending  user- 
supplied  text  to  archive  files,  and  the  like. 


1 1 
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APPENDIX  A 
CLUSTERING  ALGORITHMS 


A. 1  A  ONE-PASS  ALGORITHM 

The  -following  algorithm  works  well  when  the  data  are  grouped 
into  clusters  whose  diameter  is  less  than  the  distance  to  their 
nearest  neighboring  cluster.  Signi-ficant  violation  o-f  this 
condition  will  make  the  algorithm  highly  sensitive  to  the  arbi¬ 
trary  choices  imposed  in  the  ordering  of  the  data. 

In  the  tollowing  description,  the  data  points  X,  Y,  etc., 
maxy  be  taken  to  be  scalar  voltage  samples,  in  the  case  where  we 
are  looking  for  clusters  in  the  amplitude  data;  or  they  may  be 
taken  to  be  the  vectors  of  dimension  N  consisting  of  the  first 
N  samples  including  and  following  a  threshold  crossing,  in  the 
case  where  we  are  looking  for  clusters  in  the  pulse  waveform 
types.  In  either  case,  the  notation  !X-Y!  means  the  Euclidean 
distance  between  the  points,  whether  they  are  in  one  dimension  or 
in  N  dimensions. 

Let  (XI , X2, . . . , Xn)  be  a  sequence  of  n  data  points  (or 
vectors),  and  define  Zi  to  be  the  i-th  cluster  center.  The 
algorithm  finds  the  set  fZi>  of  cluster  centers.  First,  it  is 
required  to  obtain  a  cluster  diameter,  D,  and  we  assume  here  that 
we  can  obtain  it  by  experimentation  (to  see  which  values  produce 
the  most  reasonable  clusterings)  or  by  a-priori  (knowledge  of  the 
variance  in  amplitudes  from  a  single  emitter. 

Having  established  D  we  then  define  Zl  =  XI.  Then,  for 
i  =  2  to  n  ,  compute  the  distance  from  Xi  to  each  of  the  cluster 
centers,  Zj.  If  the  distance  is  greater  than  D  for  each  j, 
then  add  Xi  to  the  set  of  cluster  centers.  Otherwise,  assign 
Xi  to  the  first  cluster  for  which  !Xi-Zjl  <  D. 

If  the  clusters  are  sufficiently  well-defined  that  this  is  a 
reasonable  algorithm  to  use,  then  we  recommend  a  post-clustering 
step  to  redefine  the  cluster  centers  (which  will  be  used  later  for 
the  realtime  amplitude  di scr i mi nat i on )  to  be  the  centroids  of  the 
individual  clusters,  found  by  vector  or  scalar  averaging. 


A. 2  The  MAXIMIN  ALGORITHM 

The  Maximin  clustering  algorithm,  like  the  one  described 
above,  is  a  heuristic  procedure  but  in  this  case  multiple  passes 
through  the  data  are  required.  The  difference  is  that  the  Maximin 
looks  for  clusters  that  are  farthest  apart  first,  and  instead  of 
having  to  know  an  explicit  feature  of  the  clusters  in  advance 
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(i.e.,  the  diameter,  D) ,  we  have  to  ex pet i men t  with  a  more 
nebulous  parameter,  F,  which  is  a  number  between  O  and  1. 

Start  with  F  =  0.5  -for  now. 

As  with  the  one-pass  algorithm,  we  begin  by  selecting  the 
•first  datum  (scalar  or  vector)  XI  to  be  the  first  cluster  center, 
21.  For  the  second  cluster  center,  22,  we  find  the  data  point 
that  is  farthest  from  21.  (If  the  distance  between  Z1  and  22  is 
sufficiently  small,  we  can  declare  that  there  is  only  one  cluster 
and  quit.)  Let  A1  be  the  distance  ! 22-21!. 

Suppose  we  now  have  a  set  C21,...,2mj'  of  cluster  centers,  a 
number  A(m-l)  which  is  the  average  of  the  previous  maximum 
distances,  and  let  (Yl,...,Ynj'  be  the  set  of  data  points  that  have 
NOT  been  assigned  to  clusters  yet.  For  each  j  =  l...m,  compute 
the  distances  Di j  =  !2j-Yi!,  i  =  l...n  ,  and  save  the  MINIMUM  of 
these,  say,  Dj’.  Then  find  the  MAXIMUM  of  the  fDj'  :  j=l,...,m>. 
Call  it  D’ .  If  D'  is  greater  than  F#A(m-l)  then  declare  the 
sample  corresponding  to  D''  to  be  a  new  cluster  center,  2(m'<-l), 
and  compute  the  new  average  maximum  distance  with  D'  included. 
Otherwise,  terminate  the  algorithm. 


1  4 


MARTINGALE  RESEARCH  CORPORATION 
NTSU  TECHNICAL  REPORT  8/01/86 


APPENDIX  B 


PEAK-FINDING  ALGORITHMS 


B. 1  THE  GRADIENT-SIGN -CHANGE  ALGORITHM 

It  is  well-known  Irom  elementary  calculus  that  an  e>:tremum  o-f 
a  differentiable  function  occurs  wherever  the  siqn  of  the 
derivative?  fgraclient)  changes  from  positive  to  negative,  or  vice 
versa.  This  fact  has  been  used  by  numerous  authors  (cf.  C4D),  and 
it  provides  a  fast  and  easy  computational  method  whenever  the 
signal  exhibits  wel 1 -separated  peaks  that  are  well  above  the 
noise.  Those  conditions  seem  to  apply  in  the  present  situation. 

The  algorithm  is  applied  whenever  the  thresholding  detector 
has  declared  that  the  signsil  has  crossed  the  threshold,  and  it 
continues  until  an  end  condition  is  satisfied.  Let  XI  be  the 
datum  whose  absolute  value  has  just  exceeded  the  threshold,  and 
let  n  be  the  least  number  of  samples  that  are  ever  needed  at  the 
present  sample  rate  to  cover  any  spike  waveform  in  the  data.  For 
each  j  =  1  ...  n  we  check  that  the  sign  of  <Xj-X(j-l))  is 
different  from  the  sign  of  (X(j+1)-Xj).  Zero  is  included  as  a 
possible  third  "sign".  If  the  sign  has  changed,  then  an  extremum 
of  height  (or  depth,  if  negative)  Xj  is  declared  to  have 
occurred  at  the  index  (time)  j,  and  the  current  index  is  given  an 
increment  (in  ADDITION  to  the  normal  loop?  increment)  so  that  in 
case  (X(j+1)-Xj)  was  zero  the  next  point  will  not  also  be 
declared  as  a  peak. 

The  algorithm  continues  until  a  predetermined  number  of  peaks 
have  been  found,  or  until  the  n-th  datum  following  the  threshold 
crossing  has  been  tested,  whichever  occurs  first.  Processing  the 
n-th  datum  without  finding  a  peak  must  be  reported  as  an  error. 

The  manner  of  detecting  that  the  two  differences  have  changed 
sign  depejnds  on  the  number  of  CPU  clock  cycles  that  are  required 
for  a  Tiultiply  instruction.  If  the  product  of  the  adjacent 
differences  is  less  than  or  equal  to  zero,  then  the  sign  has 
changed.  However,  if  a  multipily  is  too  costly,  then  the  logical 
decisions  can  be  streamlined  by  keeping  track  of  whether  the 
threshold  crossing  was  a  negative  crossing  or  a  positive  crossing 
and  using  one  of  two  sequences  of  logical  tests,  one  being 
optimized  for  ascending  data,  the  other  for  descending  data,  with 
a  switch  being  made  after  each  peak  is  found. 

The  reason  that  one  might  want  to  find  more  than  one  peak  in 
the  waveform  is  that  it  provides  an  additional  parameter  for 
ampliti.(de  di  scr  i  mi  nati  on  that  might  be  Lised  successfully  to  avoid 
tiaving  to  classify  the  peak  with  a  Euclidean  metric  in  20  to  50 
d  i  men si  on s . 
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